Gravitational- Wave Extraction from Neutron-Star Oscillations: 
comparing linear and nonlinear techniques. 
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The main aim of this study is the comparison of gravitational waveforms obtained from numerical 
simulations which employ different numerical evolution approaches and different wave-extraction 
techniques. For this purpose, we evolve an oscillating, nonrotating, polytropic neutron-star model 
with two different approaches: a full nonlinear relativistic simulation (in three dimensions) and 
a linear simulation based on perturbation theory. The extraction of the gravitational-wave signal 
is performed via three methods: the gauge-invariant curvature-perturbation theory based on the 
Newman- Penrose scalar %pi\ the gauge-invariant Regge-Wheeler-Zerilli-Moncrief metric-perturbation 
theory of a Schwarzschild space-time; some generalization of the quadrupole emission formula. 
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I. INTRODUCTION 

The computation of the gravitational-wave emission 
from compact sources like supernova explosions, neutron- 
star oscillations and the inspiral and merger of two 
compact objects (like neutron stars or black holes) is 
one of the most lively subjects of current research in 
gravitational-wave astrophysics. This goal may be pur- 
sued using different numerical approaches. That is, (i) 
solving the full set of coupled Einstein and matter equa- 
tions; (ii) solving the linearized Einstein and matter equa- 
tions around a fixed background, when such an approx- 
imation is valid. In the latter case, with the additional 
condition of spherical symmetry, the formalism we em- 
ploy is based on a multipolar expansion and the compu- 
tation of the gravitational waves directly follows from the 
knowledge of the pcrturbative metric multipoles k? m , xtm 
and ipi m . On the other hand, extracting gravitational 
waveforms from a space-time computed numerically in a 
given coordinate system is a highly nontrivial problem 
that has been addressed in various ways in the litera- 
ture. In general, two routes have proven successful: (i) 
the gauge-invariant curvature-perturbation theory based 
on the Newman-Penrose [1] scalar ^4, and (ii) the Regge 
and Wheeler [2] , Zerilli [3] theory of metric-perturbations 
of a Schwarzschild space-time, recast in a gauge- invariant 
framework following the work of Moncrief [4] . 

The aim of our study is the computation of the grav- 
itational waveforms emitted by the very controlled sys- 
tem constituted by a nonrotating polytropic relativistic 
star that oscillates nonisotropically around its spherically 
symmetric equilibrium configuration because of an ax- 
isymmetric perturbation. Our aim is to follow two (com- 
plementary) calculation procedures. On one hand, wc 
perform a full 3+1 numerical simulation of the system, 



i.e. we compute a numerical solution of the Einstein equa- 
tions without approximations except those of the numeri- 
cal method itself. Because of its generality, this approach 
allows us to analyze different physical regimes, in partic- 
ular, the case in which the "perturbation" is not small 
and nonlinear effects can play a relevant role with impor- 
tant consequences on the waveforms. On the other hand, 
we follow a perturbative approach based on the assump- 
tion that the perturbation is "small" . If this is the case, 
one can (i) expand the metric around a fixed background 
(i.e. the Tolman-Oppenhcimcr-Volkoff solution), (ii) re- 
tain only the linear term of this expansion and (iii) solve 
the linearized Einstein equations. In addition, since the 
star is nonrotating, one can factorize the angular depen- 
dence by means of a spherical-harmonic decomposition 
of the metric and matter fields, and, thus, only a 1+1 
system of partial differential equations must be solved. 

The present work has much in common with Refs. [5, 
6] , where a comparison of different extraction techniques 
has been performed. Following the same inspiration of 
Ref. [6], we exploit perturbative computations to ob- 
tain "exact" waveforms to compare with the numcrical- 
relativity-generated ones. As done in Ref. [5], we use 
an oscillating neutron star as a test-bed system, but we 
consider a wider range of possible wave-extraction tech- 
niques. Since there is a copious literature dealing with 
the problem of gravitational-wave extraction in numer- 
ical relativity, we prefer not to mention here the main 
bibliographic references, but rather to address the reader 
to the references in Refs. [5, 6] and to the citations in the 
following text. 

The article is organized as follows. In Sec. II we 
describe the numerical time-evolution methods and the 
gravitational-wave extraction techniques adopted. In 
Sec. Ill we introduce our choice of initial data and Sec. IV 
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is devoted to the presentation of our results. Conclu- 
sions that can be drawn from our results are discussed in 
Sec. V. 

Standard dimcnsionless units c = G = M@ = 1 and a 
spacelike signature (— , +, +, +) are used. Greek indices 
are taken to run from to 3, Latin indices from 1 to 3 
and we adopt the standard convention for the summation 
over repeated indices. 



II. THE PHYSICAL SYSTEM AND ITS 
NUMERICAL EVOLUTION 

In this section we present the main elements of the 
two evolutionary approaches and discuss the three wave- 
extraction techniques mentioned in the introduction. In 
our investigation we deal with the full set of Einstein 
equations 



(1) 



coupled to a perfect-fluid matter, with stress-energy ten- 
sor 



T» v = p ! + £+'- \u»u v +pg< iV , 



(2) 



where u M is the fluid 4- velocity, p is the fluid pressure, 
e is the specific internal energy and p is the rest-mass 
density, so that e = p(l + e) is the energy density in 
the rest frame of the fluid and H = p(\ + e) + p is the 
relativistic specific enthalpy. The Einstein equations for 
the space-time must be supplemented by the relativistic 
hydrodynamics equations, namely, the conservation law 
for the energy-momentum tensor V M T'"' = 0, the con- 
servation law for the baryon number V A1 (/5M A ') = 0, and 
an equation of state (EOS) of the type p = p(p, e). For 
the purpose of this work, we restrict our attention to the 
polytropic (isoentropic) equation of state: 



p = Kp L 
K 



,r-i 



(3) 



r-i' 



with parameters K = 100 and T = 2. 



perturbations of spherically symmetric space-times devel- 
oped in Refs. [12-16]. We work explicitly in the Regge- 
Wheeler gauge. In this case, the full set of perturbation 
equations that we use is equivalent to that of Refs [17, 18]. 

The focus of this work is on even-parity perturba- 
tions only 1 . Let us recall that Ref. [9] showed how the 
even-parity perturbation problem can be set up, and sta- 
bly solved, using a constrained formulation of the per- 
turbation equations. These equations, as well as their 
numerical solution, have been discussed several times 
in the literature [8, 9, 11]. Notably, common practice 
is that (i) one elliptic equation, the Hamiltonian con- 
straint, namely Eq. (7) of Ref. [11], is solved to ob- 
tain the perturbed conformal factor, kg m \ (ii) one hy- 
perbolic equation, namely Eq. (6) of Ref. [11], is used 
(only inside the star) to evolve the matter variable Hg m 
(i.e. the perturbation of the relativistic enthalpy); (iii) 
another hyperbolic equation, namely Eq. (5) of Ref. [11], 
permits to obtain the nondiagonal, gauge-invariant met- 
ric degree of freedom (the one actually associated with 
gravitational radiation), xim- After specification of ini- 
tial data, the hyperbolic equations are solved with stan- 
dard, second-order-convergent-in-time-and-space, finite- 
differencing algorithms (e.g. leapfrog or Lax-Wendroff ). 
Consistently, the elliptic equation is discretized at second 
order in space and reduced to a tridiagonal linear system, 
which is then solved by inversion. For any given multi- 
pole, (£, m), one solves the system of equations to obtain 
Xim and kg m as functions of time. Outside the star, one 
finally computes the Zerilli-Moncrief function as 



(e) 



2r(r - 2M) 
A[(A- 2)r + 6M] 



Xtr, 



rd r kir, 



rA + 2M 
2(r - 2M) 



tilt 



(4) 



where M is the stellar mass and A = £(£+1). This func- 
tion is directly connected to the h + and h x gravitational- 
wave polarization amplitudes [sec Eq. (40) below] and it 
can be extracted from general-relativistic 3D codes; for 
this reason it will be the main object of our interest in 
the forthcoming discussion. Note that Eq. (4) also de- 
fines our normalization conventions and notation, that 
agree with those of Ref. [19]. 



A. PerBACCo: a general-relativistic ID linear code 



The PerBACCo (PerturBAtive Constrained Code) 
general-relativistic linear code that we employ in this 
work is a development of the one introduced in Refs. [7, 8] 
and recently used in many studies [9-11]. This code is 
1+1-dimcnsional and evolves, in the time domain, non- 
spherical, matter and metric linear perturbations of a 
spherical star. The equations that are solved are ob- 
tained, after a multipolar decomposition of the linearized 
Einstein equations, as the static-background case in the 
gauge- invariant and coordinate-independent formalism of 



B. Cactus-Carpet-CCATIE-Whisky: 



1 The metric perturbations of a spherically symmetric background 
space-time are divided in two classes, which are decoupled: the 
even-parity perturbation (also called electric because it is gener- 
ated by the time variation of the mass multipole moments of the 
source), which transform as (— 1)* under a parity transformation, 
and the odd-parity perturbation (also called magnetic because it 
is generated by the current multipole moments), which transform 
as 
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a general-relativistic 3D nonlinear code 

We evolve a conformal-traceless "3 + 1" formulation of 
the Einstein equations [20-23], in which the space-time 
is decomposed into three-dimensional spacclikc slices, de- 
scribed by a metric 7^, its embedding in the full space- 
time, specified by the extrinsic curvature K+j, and the 
gauge functions a (lapse) and (3 l (shift), that specify a 
coordinate frame (see Sec. II B 1 for details on how we 
treat gauges and Ref. [24] for a general description of 
the 3+1 split). The particular system which we evolve 
transforms the standard ADM variables as follows. The 
three-metric 7^ is conformally transformed via 



1 

12 



In dct 7jj , 



(5) 



and the conformal factor <f> is evolved as an indepen- 
dent variable, whereas 7^ is subject to the constraint 
det7ij = 1. The extrinsic curvature is subjected to the 
same conformal transformation and its trace tr is 
evolved as an independent variable. That is, in place 
of Kij we evolve 



K = tx K tJ = //'*' A, . , Aij = e-^(K l3 - -^K), 
with tr Aij = 0. Finally, new evolution variables 



(6) 



(7) 



are introduced, defined in terms of the Christoffel sym- 
bols of the conformal three- metric. 

The Einstein equations specify a well-known set of evo- 
lution equations for the listed variables and are given by 



-p) Ho 



~2aA, 



(8) 
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aK, 



(9) 



(d t - C ) Atj = e-^[-DiD ja + a(i? y - - 8tt%)] tf 

+ a{KA tj 
(d t -C P )K= -D l D ia 



a(KA l} - 2A ik A k j), (10) 



a 



(11) 



'),\" = j^fydkF + -fidjdkp* + fill' - vn r f 

+ ^FdjF - 2A ij dja + 2a(T i jk A jk + GA^dj 



(12) 



where i?^ is the three-dimensional Ricci tensor. Di the 
covariant derivative associated with the three- metric 7^ . 



"TF" indicates the trace-free part of tensor objects and 
p , DM , S j , and Sij are the matter source terms defined as 



Padm = n a n f3 T a ' 3 , 
Si = - lia n T afi , 



(13) 



where n a = (—a, 0,0,0) is the future-pointing four- 
vector orthonormal to the spacelike hypersurface and 
T a/3 is the stress-energy tensor for a perfect fluid [cf. 
Eq. 2] . The Einstein equations also lead to a set of phys- 
ical constraint equations that are satisfied within each 
spacelike slice: 

n = r {3) 



0, 



M* = Dj(K ij - j i] K) - 8ttS 1 = 0, 



(14) 
(15) 



which are usually referred to as Hamiltonian and momen- 
tum constraints. Here = i?y7 u is the Ricci scalar 
on a three-dimensional time-slice. Our specific choice of 
evolution variables introduces five additional constraints, 



det7y = 1, 



jk- 



(16) 
(17) 
(18) 



Our code actively enforces the algebraic constraints (16) 
and (17). The remaining constraints, Ti 1 and (18), 
are not actively enforced and can be used as monitors of 
the accuracy of our numerical solution. See Ref. [25] for 
a more comprehensive discussion of the above formalism. 



1. Gauges 

We specify the gauge in terms of the standard ADM 
lapse function a, and shift vector /3* [26]. We evolve the 
lapse according to the "1 + log" slicing condition [27]: 



d t a - f5 l d t a = -2a(K - K ), 



(19) 



where K is the initial value of the trace of the extrinsic 
curvature and equals zero for the maximally sliced initial 
data we consider here. The shift is evolved using the 
hyperbolic T-driver condition [25], 



a, ir - vnjr 



-aB* , (20) 

a t r - pep - r)B i , (21) 



where 77 is a parameter which acts as a damping coef- 
ficient. The advection terms on the right-hand sides of 
these equations have been suggested in Rcfs. [28-30]. 

All the equations discussed above are solved using 
the CCATIE code, a three-dimensional finite-differencing 
code based on the Cactus Computational Toolkit [31]. 
A detailed presentation of the code and of its conver- 
gence properties have been recently presented in Ref. [32]. 
Mesh refinement is achieved through the Carpet code 
[33]. 
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2. Evolution system for the matter 

We solve the general-relativistic hydrodynamics equa- 
tions with the Whisky code [34-38]. An important fea- 
ture of the Whisky code is the implementation of a flux- 
conservative formulation of the hydrodynamics equa- 
tions [39-41], in which the set of conservation equations 
for the stress-energy tensor T Ml/ and for the matter cur- 
rent density J M = pu^ ', namely 

V M T^ = 0, V M J M = 0, (22) 

is written in a hyperbolic, first-order and flux- 
conservative form of the type 

d f q + dT«(q) = s(q) , (23) 

where f^(q) and s(q) are the flux vectors and source 
terms, respectively [42]. Note that the right-hand side 
(the source terms) does not depend on derivatives of the 
stress-energy tensor. Furthermore, while the system (23) 
is not strictly hyperbolic, strong hypcrbolicity is recov- 
ered in a flat space-time, where s(q) = 0. 

As shown by Ref. [40], in order to write system (22) 
in the form of system (23), the primitive hydrodynami- 
cal variables (i.e. the rest- mass density p, the pressure p 
measured in the rest-frame of the fluid, the fluid three- 
velocity v 1 measured by a local zero-angular momentum 
observer, the specific internal energy e and the Lorentz 
factor W) are mapped to the so-called conserved vari- 



ables q = (D, S 1 ,t) via the relations 

D ee ^jWp , (24) 

S l = ^ P HW 2 v l , (25) 

r = ^(pHW 2 -p) -D . (26) 



Note that, in the case of a general EOS of the type 
V — p{Pi e ) only five of the seven primitive variables are 
independent. Furthermore, if one adopts - as we do in 
the present work - a simpler isoentropic EOS of the type 
p = p(p) where also the specific energy (e) is fully deter- 
mined by the rest-mass density (p), there is even one less 
independent variable. Namely Eq. (26) becomes redun- 
dant and needs not be solved. No fundamental changes 
need being applied to the code, except that a simpler con- 
version scheme from conservative variables to primitive 
variables can be adopted [43, 44]. 

In this approach, all variables q are represented on 
the numerical grid by cell- integral averages. The func- 
tions that the q represent are then reconstructed within 
each cell, usually by piecewise polynomials, in a way 
that preserves conservation of the variables q [45] . This 
operation produces two values at each cell boundary, 
which are then used as initial data for the local Rie- 
mann problems, whose (approximate) solution gives the 
fluxes through the cell boundaries. A method-of-lines ap- 
proach [45], which reduces the partial differential equa- 
tions (23) to a set of ordinary differential equations 
that can be evolved using standard numerical methods, 



such as Runge-Kutta or the iterative Cranck-Nicholson 
schemes [46, 47], is used to update the equations in 
time (see Ref. [48] for further details). The Whisky 
code implements several reconstruction methods, such as 
total-variation-diminishing (TVD) methods, essentially- 
non-oscillatory (ENO) methods [49] and the piecewise 
parabolic method (PPM) [50] . Also, a variety of approx- 
imate Riemann solvers can be used, starting from the 
Harten-Lax-van Leer-Einfeldt (HLLE) solver [51], over 
to the Roe solver [52] and the Marquina flux formula [53] 
(see Ref. [38, 48] for a more detailed discussion). 

In this work we always use a global second-order ac- 
curate scheme, where time evolution is performed using 
the Iterative Cranck-Nicholson scheme with three sub- 
steps and with a Courant-Friedrichs-Lewy factor equal 
to 0.25. We always use the PPM method (that it is 
nominally 3rd-order accurate, but in actual simulations 
usually shows at best second-order accuracy) for the re- 
construction and the Marquina formula for the approx- 
imate fluxes. The employed finite differencing for the 
space-time evolution with the CCATIE code is fourth- 
order accurate. There are no particular reasons to prefer 
these schemes with respect to others used in the litera- 
ture (like 3rd-order Runge-Kutta methods for time evo- 
lutions), however, since in this work we have focused on 
comparing gravitational- wave-extraction methods rather 
than time-evolution methods, we decided to use the old- 
fashioned iterative Cranck-Nicholson scheme. 



3. Treatment of the atmosphere 

At least mathematically, the region outside our initial 
stellar models is assumed to be perfect vacuum. Inde- 
pendently of whether this represents a physically realistic 
description of a compact star, the vacuum represents a 
singular limit of the Eqs. (24-26) and must be treated ar- 
tificially. We have here followed a standard approach in 
computational fluid-dynamics and added a tenuous "at- 
mosphere" filling the computational domain outside the 
star. 

We treat the atmosphere as a perfect fluid governed 
by the same polytropic EOS used for the bulk matter, 
but having a zero coordinate velocity. Furthermore, its 
rest-mass density is set to be several (6 in the present 
case) orders of magnitude smaller than the initial central 
rest-mass density. 

The evolution of the hydrodynamical equations in grid- 
zones where the atmosphere is present is the same as the 
one used in the bulk of the flow. Furthermore, when the 
rest mass in a grid-zone falls below the threshold set by 
the atmosphere, that grid-zone is simply not updated in 
time and the values of its rest-mass density and velocity 
are set to those of the atmosphere. 
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C. Gravitational- wave extraction in 
Cactus-Carpet-CCATIE-Whisky 

On a flat space-time, it is natural to express the wave- 
form as a multipolar expansion in spin- weighted spherical 
harmonics of spin weight s = — 2 as 

h+ - ifc x = £ E h im - 2 Y im {6, <f>). (27) 

i=2 m=-l 

The problem of gravitational-wave extraction out of a 
space-time computed numerically amounts to comput- 
ing, in a coordinate-independent way, the multipolar co- 
efficients h tm . Two routes are commonly followed in 
numerical-relativity simulations of astrophysical systems 
which do not involve matter (like binary black-hole co- 
alescence). On one hand, one focuses on Weyl "curva- 
ture" waveforms [54], by extracting from the numerical 
space-time the Newman-Pcnrose scalar ipi, which is re- 
lated to the second time derivative of (h+,h x ) (see be- 
low). The metric waveform (27) is then obtained from 
the curvature waveform via time integration. On the 
other hand, one can rely on the Reggc- Wheeler [2] and 
Zerilli [3] theory of metric perturbations of Schwarzschild 
space-time, after recasting it in its gauge-invariant form 
according to Moncrief [4]. This allows to compute the 
metric waveform directly from the numerical space-time. 
See also Refs. [19, 55, 56] for reviews and generalizations. 
Moreover, if matter is involved, it is also possible to cal- 
culate the gravitational radiation emitted by the system 
by means of some (modified) Landau-Lifshitz quadrupolc 
formula. The purpose of this section is to review the main 
elements of the three wave-extraction procedures, as an 
introduction to Sec. IV, where waveforms obtained via 
the different methods will be compared and contrasted. 

1. Wave- extraction via Newman- Penrose scalar 04 

The use of Wcyl scalars for wave-extraction pur- 
poses has become very common in numerical relativity 
and it has been successfully applied in current binary- 
black-hole (see Rcf. [57] and references therein), binary- 
neutron-star [34] and mixed-binary [58] simulations. 

Given a spatial hypersurface with timelike unit normal 
and given a spatial unit vector r M in the direction of 
the wave propagation, the standard definition of 04 is the 
following component of the Weyl curvature tensor C ail p v 

04 = -C a(I ^ffm m^, (28) 

where = l/v / 2(n A1 — r^) and m M is a complex null 
vector (such that m^fh^ = 1) that is orthogonal to r M 
and n M . This scalar can be identified with gravitational 
radiation if a suitable frame is chosen at the extraction 
radius. On a curved space-time there is considerable free- 
dom in the choice of the vectors r M and m M and differ- 
ent researchers have made different choices, which are 



all equivalent in the r — ► 00 limit (see for example [59] 
and references therein). We define an orthonormal basis 
in the three-space (e r ,£0,e^), centered on the Cartesian 
origin and oriented with poles along the z-axis. The nor- 
mal to the slice defines a timelike vector &t, from which 
we construct the null frame 

Z = -4=(et — e r ), n = -j=(e t + e r ), m = — (eg - ie^). 

(29) 

We then calculate 04 via a reformulation of Eq. (28) in 
terms of ADM variables on the slice [60], 

04 = C lj m l fh\ (30) 

where 

d, = Rij - KK l3 + KfKkj - ie^Vi-Kj*. (31) 

The gravitational-wave polarization amplitudes h+ and 
h x are related to -04 by [61] 

h + — ih x = ip£. (32) 

It is then convenient to expand 04 m spin- weighted spher- 
ical harmonics of weight s = — 2 as 

00 1 

Mt,r,9,<f>) = J2 E Tpl m (t,r)-2Y im (9,<f>), (33) 

t=2 m=-i 

so that the relation between 04™ and the metric multi- 
poles h lm becomes 

¥ m {t,r) = 0f"(t,r) . (34) 

h lm (t,r) is then the double indefinite integral of 
ipl m (t,r), which we numerically compute (after multi- 
plying both sides by r) as 

r~h em {t,r) = f dt' f dt"r^ 4 m {t",r) , (35) 
Jo Jo 

which results in 

r h em {t, r) = rh im (t, r) + Q Q + tQ 1 , (36) 

where the integration constants Qq and Q\ are explicitly 
written. They can be determined from the data them- 
selves and their physical meanings are Qo = rh (0,r) 
and Q x = r/0 m (O,r). 

This is not the end of the story yet. The equations 
discussed so far refer to a signal extracted at a finite value 
of r, while one is interested in computing ip^" at spatial 
infinity. It is imaginable that in the computed values 
of 0| m (0r) there may be an offset, dependent on the 
extraction radius; that is, 0| m a t spatial infinity should 
be written as 

r0f (i)^r0f l (i,r) + 2Q 2 (r), (37) 
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where ipl" l (t, r) is the scalar extracted at a finite radius r 
and 2Q 2 (r) is an offset function, that takes into account 
(in an additive way) the effects of the extraction at a 
finite radius. The time integration of this offset generates 
an additional term that is quadratic in time, so that the 
final result for rh em (t) is 

r h tm (t) = r h im (t, r) + Q + Qit + Q 2 (r)t 2 . (38) 

The term Q2 (r) should tend to zero when the extraction 
radius goes to infinity. We checked that this is the case for 
the results of our simulations (see Sec. IV B and Fig. 3). 

Various ways of fixing the two integration constants Qo 
and Qi have been discussed in the literature about coa- 
lescing binary black-hole systems [32, 62-64]. In partic- 
ular, in Appendix A of Ref. [64] the following procedure 
was proposed: (i) integrate the curvature waveform twice 
forward in time (starting from t = and including the 
initial burst of radiation due to the initial-data setup); 
(ii) Subtract the linear-in-time offset present in there. 
This simple procedure led to an accurate metric wave- 
form which exhibited the correct circular polarization be- 
havior. A similar line was also followed in Ref. [62] , where 
it was pointed out that in some situations (e.g. close ex- 
traction radius, higher multipolcs) one needs to subtract 
a general polynomial in t, consistently with our Eq. (38). 

2. Abrahams-Price metric wave-extraction procedure 

The wave-extraction formalism based on the perturba- 
tion theory of a Schwarzschild space-time was introduced 
by Abrahams and Price [65] and subsequently employed 
by many authors [66-69]. 

The assumption underlying this extraction method is 
that, far from the strong-field regions, the numerical 
space-time can be well approximated as the sum of a 
spherically symmetric Schwarzschild "background" 
and a nonspherical perturbation h^. Even if based 
on the gauge-invariant formulation of perturbations due 
to Moncrief [4], the standard implementation [65] of 
this approach is done by fixing a coordinate system 
(Schwarzschild coordinates) for the background. As 
usual, the spherical symmetry 2 of g° allows one to elim- 
inate the dependence on the angles (6>, </>) by expand- 
ing in (tensor) spherical harmonics, i.e. seven even- 
parity and three odd-parity multipoles. The multipolar 
expansion explicitly reads 

00 £ 

M*.r,M)=&+E E [(C") (0) + (C n ) (c) • 

1=2 m=-l 

(39) 



The metric multipoles (h 1 ™) ° ° (and their derivatives) 
can be combined together in two gauge-invariant mas- 
ter functions, the even-parity (Zerilli-Moncrief) [see 

Eq. (4) above] and the odd-parity (Regge- Wheeler) V E'^ . 
These two master functions satisfy two decoupled wave- 
like equations with a potential 3 . Finally, in a radiative 
coordinate system we have 

^=^(*a+i*2), (40) 

where N e = ^(£ + 2)^+1)^-1). 

Note that the use of Schwarzschild coordinates for the 
background metric is not at all necessary and more gen- 
eral wave-extraction frameworks exist. In particular, Sar- 
bach and Tiglio [55] and Martel and Poisson [56] have 
shown that there exists a generalized formalism for per- 
turbations that is not only gauge invariant (i.e. invari- 
ant under infinitesimal coordinate transformation), but 
also coordinate independent, in the sense that it is in- 
variant under finite coordinate transformations of the 
M 2 Lorentzian submanifold of the background. Since 
in a numerical-relativity simulation the gauge depends 
on time, one is a priori expecting that the gauge fixing 
of the background may introduce systematic errors. For 
the odd-parity case, Ref. [6] has shown that this is indeed 
the case for the particular physical setting represented by 
the scattering of a Gaussian pulse of gravitational waves 
on a Schwarzschild black hole in Kerr-Schild coordinates 
(see Ref. [70] for the even-parity case). In this work we 
present results obtained using the "standard" Moncrief 
formalism. A comprehensive discussion of results ob- 
tained via the generalized formalism will be presented 
elsewhere [71]. 

3. Landau- Lifshitz quadrupole-type formula 

In the presence of matter, it is sometimes convenient 
to extract gravitational waves using also some kind of 
(improved) Landau-Lifshitz "quadrupolc" formula. Al- 
though this formula is not gauge invariant, this route 
has been followed by many authors with different degrees 
of sophistication [5, 72-75], to give well approximated 
waveforms [5]. For the sake of completeness, let us re- 
view how this quadrupole formula came into being, as 
the first contribution in a multipolar expansion, and let 
us express it in the convenient form of h , as outlined 
above. The basic reference of the formalism is a review 
by Thorne [76] ; most of the useful formulas of this review 
have been collected by Kidder [77], who condenses and 
summarizes the gravitational- wave-generation formalism 
developed in Refs. [78, 79]. 



2 That is, the background 4-manifold M can be written as M = 
M 2 X S 2 , where M 2 is a two-dimensional Lorentzian manifold 3 The equations are just approximately satisfied on the extracted 
and S 2 is the unit two-sphere. background. 
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Following Rcf. [77], we recall that Eq. (40) can be 
derived in all generality by (i) decomposing the asymp- 
totic waveform hj^ into two sets of symmetric trace-free 
(STF) radiative multipole moments (to be related later 
to the matter multipole moment of the source in the near- 
zone) called IAl and Vl, where a capital letter for an index 



denotes a multi-index {i.e. , IAl = Ui 



s); (ii) project- 



ing the STF-decomposed hj^ along an orthonormal triad 
that corresponds to that of the spherical coordinate sys- 
tem. In the same notation of Ref. [77] , Eq. (40) reads 



■Am 



1 



V2r 



(41) 



where the mass multipole moments U lm and current mul- 
tipole moments V are related to their STF counter- 
parts by 



U 



V 



e m _ 167T (£ + !)(£ + 2) hw 
(2€ + 1)!! y 2£(£-l) L ' 

32n£ [ 



(42) 



im 



£ + 2 



(21 + 1)!! V 2£(£ + 1)(£ - ij VLy L n ' ( 43 ) 



where y L m are the STF spherical harmonics. These func- 
tions form a basis of the of the (2£ + l)-dimensional vec- 
tor space of STF ^-tensors; they are related to the scalar 
spherical harmonics by 



Y 



£m 



yi m N L , 



(44) 



where Ni is a component of the unit radial vector. The 
expanded form of the STF y L m is given in Refs. [76, 80] 
(see also Eq. ( A6a) of Rcf. [78] ) . In the post-Newtonian 
(PN) wave-generation formalism of Refs. [78, 79], one 
can relate in a systematic manner the radiative multi- 
pole moments (Ul,Vl) to a set of six STF source mo- 
ments (2~l, Jl, Wlj %L,yL, Zl), which can be computed 
from the stress-energy pseudotensor of the matter and of 
the gravitational field of the source. A set of two canon- 
ical source moments (A4l,Sl) can be computed as an 
intermediate step between the source moments and the 
radiative moment. Two of the source moments, the mass 
moments II and the current moment Jl are dominant, 
while the others only make a contribution starting at 
2.5 PN order and we neglect them here. In a first approx- 
imation (i.e. neglecting the nonlinear "tail-interactions" 
as well as higher-order nonlinear interactions), the L-th 
radiative moment is given by the £-th time derivative of 
the canonical moments as 



J L 



0(e 5 / 2 ), 
-0(e 5 / 2 ), 



(45) 
(46) 



where e ~ (v/c) 2 indicates some PN ordering parameter 
of the system. As a result, the computation of Ui m and 
Vi m is straightforward. As an example (that will be used 
in the following), let us focus on the £ = 2 moments of a 



general astrophysical system with equatorial symmetry. 
In this case, the (2, 1) moment is purely odd-parity, while 
the (2,0) and (2,2) are purely even-parity. Straightfor- 
ward application of what we have reviewed so far gives 



,21 



i /1287T 



r V 45 



{^Jxz ^-*~£y 



^ = lJ^(l xx -2X xy -J yy 



(47) 
(48) 
(49) 



In the harmonic gauge, in the case of small velocity and 
negligible internal stresses (i.e. in the Newtonian limit) 
one has = J d 3 xpXiXj and J7y = J d 3 xpe a i,iXjX a v b . 
The 1 PN corrections to the mass quadrupole have been 
computed in Ref. [81]. Recently, Ref. [82] included 
1 PN correction, using an effective 1 PN quadrupole 
momentum, in the gravitational-wave-extraction proce- 
dure from supernova core-collapse simulations. As a 
complementary approach, Ref. [5] proposed to "effec- 
tively" take into account possible general-relativistic cor- 
rections by inserting in Eqs. (47-49) the following effec- 
tive "quadrupole moment" defined in terms of the "co- 
ordinate rest-mass density" p* = a^/~fu°p, 



(50) 



This presents some very useful properties: (i) it is of sim- 
ple implementation and (ii) from the continuity equation 
dtp* + di(p*v l ) = 0, one can analytically compute the 
first time-derivative of the quadrupole moment, so that 
only one numerical time-derivative needs to be evaluated. 
The last property is extremely important, in fact, on data 
computed via a second-order accurate numerical scheme 
it is not possible to calculate noise-free third derivatives, 
which are needed for the gravitational-wave luminosity. 
The accuracy of a scheme based on Eq. (50) has been 
tested in Ref. [5] in the case of neutron-star oscillations 
and was subsequently used by various authors to estimate 
the gravitational-wave emission in other physical scenar- 
ios. See for example Refs. [35, 74, 83, 84]. In order to 
get some more insight on the accuracy of possible "gener- 
alized" standard quadrupole formulas (SQFs formulas), 
we have tried the strategy exploited in Ref. [85], namely 
to test some pragmatic modifications of the quadrupole 
formula and to check which one is closer to the actual 
gravitational waveform. In practice, we start with a sort 
of generalized "quadrupole moment" of the form 



X QX i X j 



(51) 



where now, instead of the "Matter density" , we use the 
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following generalized effective densities g: 




SQF 


Q 


= P, 


(52) 


SQF1 


Q 




(53) 


SQF2 


Q 




(54) 


SQF3 


Q 


o W 
= u p = —p. 


(55) 



a 



We do not think that any of the "quadrupole formu- 
las" obtained using these generalized quadrupole mo- 
ments should be considered better than the others. Note 
that none of them is gauge invariant and, indeed, the 
outcome will change if one is considering isotropic or 
Schwarzschild-like coordinates. These formulas were 
widely used in the literature and the main purpose of 
the comparison among Eqs. (52-55) is to give an idea of 
the kind of information that can be safely assessed using 
them. We will comment more on that in the discussion 
in the following Sec. IV E. 

III. INITIAL DATA 



TABLE I: Equilibrium properties of model AO. From left to 
right the columns report: central rest-mass density, central 
total energy density, gravitational mass, radius, compactness. 



Name 


Pc 


e c M 


R 


M/R 


AO 


1.28 x 10" 3 


1.44 x 10" 3 1.40 


9.57 


0.15 



TABLE II: From the enthalpy perturbation to the relative 
magnitude of the pressure perturbation for n = 0, £ = 2, 
and m = (see Fig. 1). The minimum pressure perturbation 
occurs at some value of r on the xy plane (6 — tt/2), while 
the maximum pressure perturbation is found at some value 
of r on the z axis (9 = 0). 



Name 


A 


min(Sp/p c ) 


max(i5p/p c ) 


AO 


0.001 


-0.00125 


0.00251 


Al 


0.01 


-0.01253 


0.02506 


A 2 


0.05 


-0.06266 


0.12533 


A3 


0.1 


-0.12533 


0.25067 



As a representative model for a neutron star, we choose 
a model described by a polytropic EOS [Eq. 3] with T = givcn fluid mode can bc givcn by sctting 
2, K = 100, central rest-mass density p c = 1.28 x 10 3 
and so with rest mass M ~ 1.4. This model has been 
widely used in the literature and it is known as model 

AO in Ref. [86]. Some of its equilibrium properties are "-to = A sin 

listed in Table I. 



2R 



(57) 



A. Fluid-perturbation setup 

In both the linear and nonlinear codes, setting up 
the initial data amounts to (i) solving the Tolman- 
Oppenheimer-Volkov (TOV) equations to construct the 
equilibrium configuration; (ii) fixing an axisymmetric 
pressure perturbation; and (iii) solving the linearized 
constraints for the metric perturbations. We rewrite the 
perturbative equations in terms of enthalpy perturba- 
tions because it is more convenient. 

We set up the initial pressure perturbation as an ax- 
isymmetric multipole: 

6p(r, 9) = (p + e)H m (r)Y m (9) , (56) 

and then one is free to specify a profile for the relativistic 
enthalpy Hgo(r). Actually we limit our study to £ = 2 
(quadrupole) perturbations. Since we aim at a compari- 
son between waveforms and not at exploring the physics 
of the process of neutron-star oscillations, for our purpose 
the best system is represented by a star oscillating pre- 
cisely at one frequency, i. e. such that Hi q (r) corresponds 
to an eigenfunction of the star. We set a profile of Hta(r) 
that excites, mostly, the / mode of the star (with a small 
contribution from the first overtone). In general, as sug- 
gested in Ref. [8], an "approximate eigenfunction" for a 



where n is an integer controlling the number of nodes of 
Hio(r), A is the amplitude of the perturbation and R is 
the radius of the star in Schwarzschild coordinates. The 
case n = has no nodes (i.e. no zeros) for < r < R; 
as a result, the / mode is predominantly triggered (as 
in Ref. [8]) and the p-mode contribution is negligible. If 
n = l, the / mode is still dominant, but a nonncgligiblc 
contribution of the p\ mode is present. If n = 2, in 
addition to the fundamental and the first pressure modes 
also the pi mode is clearly present in the signal. For 
higher values of n more and more overtones are excited. 

In the following, we use the same setup, Eq. (57), 
to provide initial data in both the linear and nonlinear 
codes. Correspondingly, the computation of dp is needed 
to get a handle on the magnitude of the deviation from 
sphericity. The best indicator is given by the ratio Sp/p c , 
where p c is the central pressure of the star. Fig. 1 dis- 
plays the profile of Sp/p c , at the pole and at the equator, 
as a function of the Schwarzschild radial coordinate r 
for A = AO = 0.001. For simplicity, we consider only 
71 = perturbations, with four values for the amplitude, 
namely A = [0.001, 0.01, 0.05, 0.1], in order to see, in 
the 3D code, how the transition from linear to nonlinear 
regime occurs. Maxima and minima of the initial pres- 
sure perturbation for the different values of the initial 
perturbation amplitude A can be found in Table II. 
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FIG. 1: Profile of the relative pressure perturbation 5p/p c 
[computed on the z-axis (8 = 0) and on the xy-pla,ne (0 — 
7r/2)] obtained from Eqs. (56-57) with n — and A = AO. 



B. Metric-perturbation setup: the ID linear code 

Let us turn now to discuss the implementation of 
Eq. (57) in the two codes and the corresponding treat- 
ment of the related initial metric perturbation. As dis- 
cussed in Rcfs. [14, 15], the even-parity metric perturba- 
tion of a general (nonstatic) spherically symmetric space- 
time is described by 3 degrees of freedom (ki m , Xim and 
iptm) that are the solution of three coupled partial differ- 
ential equations. On the static TOV background, only 
ki m and Xim are independent degrees of freedom of the 
gravitational field, and their evolution equations are de- 
coupled from that ofip£ m , which can be obtained at every 
time-step once xim and ki m are known. The recovery of 
ipe m from Xim and ki m , that is needed in the 3D case, 
will be explicitly discussed in Sec. IIIC below. By con- 
trast, for the ID implementation one only needs to spec- 
ify initial data for Xtm, k(, m , and their time derivatives. 
This is accomplished by solving the constraints under a 
number of assumptions related to the physics that we 
want to investigate. First of all, we consider only ax- 
isymmetric perturbations (m = 0) and we restrict our- 
selves to the dominant quadrupole mode (I = 2). Then, 
since in this work we are not interested in w-mode excita- 
tion, we impose the conformally flat condition (x2Q = 0) 
(see Refs. [9, 11] for details). With These hypotheses, 
we solve the Hamiltonian constraint, namely Eq.(7) of 
Rcf. [11], for ^20- This is done on a grid r g [0, r max ], 
with r max 3> R and with boundary condition fe2o = at 
r = and at r = r max . We impose /c20 = X20 = for sim- 
plicity, but we are aware that this is inconsistent with the 
condition that H20 =/= initially and thus the momentum 
constraints should also be solved. However, since the 
effect is a small initial transient in the waveforms that 
quickly washes out before the quasiharmonic oscillation 



triggered by the perturbation H20 sets in, we have de- 
cided to maintain the initial-data setup simple. Figure 2 
synthesizes the information about the initial data. The 
top panel shows (as a solid line) the profile of ^20 (versus 
Schwarzschild radius) corresponding to the perturbation 
A = AO of Table II; the bottom panel shows (as a solid 
line) the initial profile of the Zcrilli-Moncricf function 

^20 outside the star. 



C. Metric-perturbation setup: the 3D nonlinear 
code 

In the 3D code we setup the same kind of initial con- 
dition as in the ID code, but the procedure is more com- 
plicated as one needs to reconstruct the full 3D metric 
on the Cartesian grid. In addition, the main difference 
with respect to the ID case is that the perturbative con- 
straints are expressed using a radial isotropic coordinate 
f instead of the Schwarzschild-likc radial coordinate r. 
This is done because f is naturally connected to the 
Cartesian coordinates in which the code is expressed, 
i.e. f = \J x 1 + y 2 + z 2 . The initialization of the met- 
ric in the 3D case has to follow four main steps: (i) The 
perturbative constraints are solved, (ii) The multipolar 
metric components are added to the unperturbed back- 
ground TOV metric; (iii) The resulting metric is written 
in Cartesian coordinates; (iv) It is interpolated on the 
Cartesian grid. 

Let us then recall some useful formulas. At the back- 
ground level, the TOV metric in isotropic coordinates 
reads 



e 2a dt 2 + e 2b (df 2 +f 2 dn 2 ) 



(58) 



where df2 = d6 2 + sin 2 dip 2 . The relations between the 
Schwarzschild and isotropic radial coordinates in the ex- 
terior are given by 



r = r 1 + 



M 
2F 



= i (Vr 2 - 2Mr + r- m) , 



and in the interior by 

,26(r) 



r = re 



r = Cr exp 



1 - sj\ - 2m{x)/x 
xy/l — 2m(x)/x 



dx 



(59) 
(60) 

(61) 
(62) 



where 



C= (VR 2 - 2MB. + R — M 
2R V 



x exp 



R , 1 - Jl - 2m(x)/a 

dx 

xyl — 2m(x)/x 



(63) 
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FIG. 2: (color online) Initial-data setup in the ID and 3D 
codes. Top panel: Profiles of the feo multipole at t = 
versus the Schwarzschild radial coordinate r, obtained from 
the solution of the perturbative Hamiltonian constraint with 
A = AO and n = 0. Bottom panel: Profiles of $20 a ' ' - 
versus the Schwarzschild radial coordinate r, for different 
values of the initial perturbation A. Both panels compare the 
results of the computations of the ID and 3D codes. 



In terms of the isotropic radius, the perturbative Hamil- 
tonian constraint explicitly reads 



,.2J> 



„2b 



1+ [ 

A 



2m(r) 



1/2' 



tin 



2m(r) 



1/2 

£m 
A , r 



A + 2 

2r 2 



- 87re ) X lm = - 



8tt(p + e)H 



tin 



C: 



(64) 



where r = r(f) according to Eq. (59) and (61). 

After solving the TOV equations, we choose a pro- 
file for Hi m , impose the conformal- flatness condition 



Xlm = 0, and solve Eq. (64) for kg m . As we discussed 
in the previous section, one also needs to impose on ipim 
some condition that can be regarded as an initial gauge 
condition. Then (xtm, ke m , ipim) must be inserted in the 
explicit expression of the (even-parity) metric perturba- 
tion 



(65) 



{(Xe m + k em )e 2a dt 2 - 2ip em e a+b dtdf 

+ e 2b [(Xlm + k (m )df 2 + f 2 k lm dQ] }Y im 

In the absence of azimuthal and tangential velocity 
perturbations, in Schwarzschild coordinates and in the 
Rcggc- Wheeler gauge, from the momentum-constraint 
equation, namely Eq. (95) of Rcf. [15], one obtains 



4>t 



2 [m(r) + 4:irr 3 p(r)] 



2m(r) 



which, once solved, gives 



ipim {r) = C exp — / dx 



2 (m(x) + 4irx 3 p(x)) 
x — 2m(x) 



(66) 



(67) 



The requirement ipi m — > for large r, like for ki m , im- 
plies C = 0; therefore, the metric perturbation is given by 
Eq. (65) with ipg m = xim = 0. The full metric in isotropic 
coordinates is obtained as ds 2 — ds^ + 8s 2 m . This met- 
ric is transformed to Cartesian coordinates and then it 
is linearly interpolated onto the Cartesian grid used to 
solve the coupled Einstein-matter equations numerically. 
To ensure a correct implementation of the boundary con- 
ditions (i.e. ki m — > when f — > oo), the isotropic radial 
grid used to solve Eq. (64) is much larger (f ~ 3000) 
than the corresponding Cartesian grid (f ~ 208) and the 
spacing is much smaller. 

One proceeds similarly for the matter perturbation: 
From a given profile for H^ m {r), the pressure perturba- 
tion 8p{r, 9) is computed and from this the total pressure 
is given by p + Sp. This is interpolated on the Cartesian 
grid to finally obtain the vector of the conserved hydro- 
dynamics variables (D,S z ,t). 

The consistency of the initial-data setup procedure 
in both the PerBACCo ID linear code and in the 
Cactus-Carpet-CCATIE-Whisky 3D nonlinear code is 
highlighted in Fig. 2. The top panel of the figure com- 
pares the profiles of &20 m the ID case (solid line) and 
in the 3D case (dashed line) for A = AO and n = 0. The 
small differences are related to a slightly different location 
of the star surface in the two setups and to the different 
resolution of the grids. The bottom panel of Figure 2 
contrasts the Zerilli-Moncrief functions ^20 fr° m the ID 
code (solid lines) with those extracted (at t = 0) from 
the numerical 3D metric (dashed lines). For all initial 
conditions, the curves show good consistency. 



IV. RESULTS 

The presentation of our results is organized in the fol- 
lowing way. In Sec. IV A we focus first on radial oscilla- 



11 



TABLE III: Frequencies of the fundamental radial mode of 
model AO. 



n 


Pert. [Hz] 


3D [Hz] 


Diff. [%] 





1462 


1466 


0.3 


1 


3938 


3935 


0.1 


2 


5928 


5978 


0.8 



tions, that are always present due to numerical discretiza- 
tion error. Then we concentrate on nonradial stellar os- 
cillations and we compare the ID and 3D metric wave- 
forms (Sec. IV B) and curvature waveforms (Sec. IV C). 
In Sec. IV D we discuss advantages and disadvantages of 
these two wave-extraction techniques. Finally, we discuss 
the use of quadrupole-type formulas in Sec. IV E, while 
Sec. IV F is devoted to the analysis of nonlinear couplings 
between oscillation modes. 



A. Radial oscillations 

The unperturbed configuration AO has been stably 
evolved for about 20 ms. The numerical 3D grid used for 
this simulation is composed of two concentric cubic boxes 
with limits [—32, 32] and [—16, 16] in all the three Carte- 
sian directions. The boxes have resolutions A xyz = 0.5 
and 0.25 respectively; bitant symmetry, i.e. the z < 
domain is copied from the z > domain instead of being 
evolved, was imposed as a boundary condition in order 
to save computational time. 

The truncation errors of the numerical scheme trigger 
(physical) radial oscillations of (mainly) the fundamen- 
tal mode F and the first overtones. We have checked 
that these frequencies agree with those computed evolv- 
ing the radial pulsation equation with the perturbative 
code. This comparison is shown in Table III. We note in 
passing that our numbers are in perfect agreement with 
those of Table I of Ref. [87]. 

As a further check, the entire sequence of uniformly ro- 
tating models with mass M = 1.4 and nonrotating limit 
AO has been evolved. Simulations were done with a cubic 
grid with limits [—32, 32] in each direction, and uniformly 
spaced with grid spacing A xyz ~ 0.5. As before, we have 
imposed bitant symmetry. The sequence of initial models 
has been computed by means of the version of the RNS 
code [88] implemented in Whisky. For the equilibrium 
properties of the models, see Ref. [89] . 

The fluid modes of this sequence were previously inves- 
tigated in different works, using various approaches [86, 
89, 90]. With our general-relativistic 3D simulations we 
are able to study the effect of rotation on the radial mode 
and compare the results with those obtained via approx- 
imated approaches. Our results are summarized in Ta- 
ble IV. We have found that the frequencies computed by 
Dimmelmeier et al. [89] in the conformally flat approxi- 
mation are consistent with ours (the difference is of the 



TABLE IV: Frequencies of the fundamental radial mode of 
models in the sequence AU of uniformly rotating polytropic 
stars of Refs. [86] and [89]. The frequency of model AU0 (AO) 
has also been computed in Ref. [87] (1450 Hz) and in this 
work (1462 Hz), where a finer grid was used (see Table HI). 
The data in the column marked as "CF" refer to Table III of 
Ref. [89]. The data in the column marked as "Cowling" refer 
to Table II of Ref. [86]. 



MODEL 


F [Hz] 


F(CF) [Hz] 


F(Cowling) [Hz] 


AU0 


1444 


1458 


2706 


AU1 


1369 


1398 


2526 


AU2 


1329 


1345 


2403 


AU3 


1265 


1283 


2277 


AU4 


1166 


1196 


2141 


AU5 


1093 


1107 


1960 



order of few percents); on the other hand, the results of 
Stergioulas et al. [86], obtained in the Cowling approx- 
imation, differ of about a factor two, consistently with 
the estimates of Ref. [91]. In all cases, the frequencies 
decrease if the rotation increases and the trend is linear 
in the rotational parameter (3 = T/\W\, the ratio be- 
tween the kinetic rotational energy and the gravitational 
potential energy. 



B. Nonradial oscillations: comparing ID and 3D 
metric waveforms 

Let us now turn to the discussion of nonspherical oscil- 
lations and to the related extraction of waveforms from 
ID and 3D simulations. We consider a star perturbed 
with an £ = 2, Heo(r) profile with n = 0, according to 
the procedure outlined in Sec. III. This system is evolved 
separately with the two codes and the related gravita- 
tional waveforms are compared. 

We focus first on the discussion of the outcome of 
the ID linear code. We accurately performed very long 
simulations, whose final time is about 1 s. The extrac- 
tion radii for the Zcrilli-Moncrief function extend as far 
as f = 420 (~ 300 M). The resolution of the radial 
grid is Ar = 0.032, which corresponds to having 300 
points inside the star. Fig. 3 shows the Zerilli-Moncrief 
function Vl/^o (f° r ^ = ^1) extracted at different radii. 
It is plotted versus the observer retarded time, namely 
u — t — r* , where r* is the Regge- Wheeler tortoise coor- 
dinate r* = r + 2Mlog[r/(2M) — 1] and M is the mass 
of the star. 

The farther observers that are shown in Fig. 3 are suf- 
ficiently deep in the wave-zone that the initial offset, that 
is typically present due to the initial profile of /c2o> is small 
enough to be considered negligible. We checked the con- 
vergence of the waves with the extraction radius using 

as a reference point the maximum of This point 

can be accurately fitted, as a function of the extraction 
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radius, with 



(• 



a 
r 



(68) 



The extrapolated quantity a°° allows an estimate of the 
error related to the extraction at finite distance 



5a 



I ( r ) 



(69) 



The values of 5a for different radii are Sa ~ 0.5 for r = 
25M, Sa ~ 0.09 for r = 50M, 6a ~ 0.017 for r = 100M 
and 5a < 0.016 for r > 200M. 

The waveform can be described by two different 
phases: (i) an initial transient, of about half a 
gravitational-wave cycle, say up to u ~ 50, related to 
the setup of the initial data 4 , followed by (ii) a quasi- 
harmonic oscillatory phase, where the matter dynamics 
are described in terms of the stellar quasinormal modes. 
From the Fourier spectrum of over a time interval 
from 1 to about 30 ms (namely u £ [50, 6000]), we found 
that the signal is dominated by the / mode (at frequency 
Vf = 1581 Hz) with a much lower contribution of the first 
p mode (at frequency around v pi = 3724 Hz). The fre- 
quency of the / mode agrees with that of Rcf. [87] within 
1 — 2%. The accuracy of our linear code for frequencies 
obtained from Fourier analysis on such long time scries 
has been checked in Refs. [8, 11] and is better than 1% 
on average. We mention that the Fourier analysis of the 
matter variable Ht m permits to capture some higher over- 
tones than the p\ mode, although they are essentially not 
visible in the gravitational-wave spectrum. In a first ap- 
proximation, the waveform can thus be thought as the 
superposition of damped harmonic oscillators 

N 

*20 ~ ^2 Mk cos(27w 2 fcW + 4>2k) cxp(-a 2 fcii) , (70) 

k=0 

and we aim at determining the quantities Azu, </>2fc, ^2k 
and ct2k from a standard nonlinear least-square fit. Since 
the frequency v is also independently known from the 
Fourier analysis, it is used as feedback for the fit. In 
addition, to quantify the global differences between the 
"actual" and the "fitted" time series, we compute the (I 2 ) 
scalar product 



e(x,Y) = 



(71) 



4 In practice, the first half cycle of the waves cannot be expressed 
as a superposition of quasinormal modes and it is related to the 
initial data setup. This initial transient is related to two facts: 
(i) We use the conformally flat approximation; (ii) We assume 
kim = even if our initial configuration (a star plus a nonstatic 
perturbation) is evidently not time symmetric, since a velocity 
perturbation is present and thus also a radiative field related to 
the past evolutionary history of the star. 



TABLE V: Perturbative ID evolution with A = Al: results 
of the fit to a superposition of two fluid modes [see Eq. (70)] 
over the interval u € [0, 200000], i.e. about 1 s (c/.Fig. 4). For 
this fit, K ~ 1.6 x 10~ 6 and V ~ 3.6 x 10~ 5 (see text for 
explanations) . 



[arbitrary units] 


v [Hz] 




a [s- 1 ] 


A2o=l-3145±2 x 10~ 3 
j4 3 i=3.517tii x 10" 5 


^20=1583.73 69;* 
V2i =3706.9413 j 


-2 
-1 
-11 
11 


020 =3.7358^ 
a 2 i=0.421±* 




sou 



FIG. 3: Evolution of the Zerilli-Moncrief function, extracted 
at various isotropic radii f, versus the retarded time u, for the 
ID linear code evolution with A = Al. Note how the initial 
offset decreases with the extraction radius. 



which is bounded in the interval [0, 1], and then we look 
at the residual 1Z = 1 — 9. This residual gives us a 
relative measure of the reliability of the fit. In addition 
we use the l°° distance: 



V(X,Y) = max|Xj - Yj 



(72) 



that gives the maximum difference between the two time 
series. We will use the quantities 1Z and T> also as mea- 
sures of the global agreement between the 3D and ID 
waveforms. 

On the interval u € [50,6000], the waves can be per- 
fectly (TZ ~ 7 x W- 4 , P~6x 10" 6 ) represented by a 
one-mode expansion, N = 1, as the waveform is dom- 
inated by /-mode oscillations. The frequency we ob- 
tain, 1^20 = 1580.79 ± 0.01 Hz, is perfectly consistent 
with that obtained via Fourier analysis; for the damp- 
ing time, we estimate CV20 = 3.984 ± 0.066 s _1 and thus 
t~2q = &2Q — 0.25 s. If we consider the entire duration 
(Is) of the signal (see the inset in Fig. 4), it is clear that a 
one-mode expansion is not sufficient to accurately repro- 
duce the waveform. The Fourier analysis of the waveform 
in two different time intervals, one for t < 0.5 s and one 
for t > 0.5, reveals that in the second part of the sig- 
nal the pi mode, which has longer damping time, clearly 
emerges and must be taken into account. We fit the en- 
tire signal with two modes, namely N = 2, with a global 
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FIG. 4: Evolution of the Zerilli-Moncrief function, extracted 
at isotropic radius f = 280 for the same evolution of Fig. 
3 versus the retarded time u expressed now in seconds (1 
dimensionless time unit is equal to 4.92549 x 10 -6 s). The 
main panel shows only the envelope of the waveform that is 
dominated by the damping time of the / mode (17 ~ 0.27 s). 
The two insets represent the full waveform at early (top) and 
late (bottom) times. The presence of the overtone is evident 
in the oscillations at late times. 



agreement of TZ ~ 2 x 1CT 6 and V ~ 4 x lfT 5 . The re- 
sults of the fit are reported in Table V. The frequencies 
are slightly larger than those computed via Fourier anal- 
ysis and via the fit procedure restricted to only one mode 
on a shorter interval. They are, however, still consistent. 
The damping times arc T20 = 0.268 s and T21 = 2.37 s, 
with errors of the order of 0.1% and 2% respectively. 

At this stage, we have clearly assessed the accuracy 
of the waveforms computed via our ID code; in the fol- 
lowing we shall consider these waveforms (extracted at 
the farthest observer) as exact for all practical purposes. 
We turn now to the discussion of the metric waveforms 
extracted from the 3D code and we compare them to 
the exact, perturbative results for different values of the 
perturbation A. 

The 3D simulations arc performed over grids with three 
refinement levels and cubic boxes with limits [—120, 120], 
[—24, 24] and [—12, 12] in each direction. The resolutions 
of each box are A xyz = 0.5, 0.25 and 0.125, respectively. 
Equations are evolved only on the first octant of the grid 
and symmetry conditions are applied. The outermost 
detector is located at isotropic-coordinate radius f = 110 
(~ 80M). 

Figure 5 is obtained with perturbation A = Al. It 
displays the Zerilli-Moncrief normalized metric wave- 
forms, extracted on coordinate spheres of radii f G 
{30, 60, 90, 110} and plotted versus the (approximate) 
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FIG. 5: Zerilli-Moncrief normalized metric waveforms shown 
versus the observer retarded time u — t — r» at different ex- 
traction radii (f = 30 to f = 110), for a 3D evolution with 
perturbation A = Al. The initial part of the waveform is 
dominated by a pulse of junk radiation. 



retarded time u = t — r* , where r* = r+2M log[r/ (2M) — 
1]. Here, r is the areal radius of the spheres of coordi- 
nate radius f and M is the Schwarzschild mass enclosed 
in f [6, 65, 69]. This figure is the 3D analogous of Fig. 3. 
The ID and 3D waveforms look qualitatively very sim- 
ilar apart from the presence of a highly-damped, high- 
frequency oscillation at early times. In Sec. IV D we will 
argue that this oscillation is essentially unphysical be- 
cause its amplitude grows linearly with the extraction 
radius f, instead of approaching an approximately con- 
stant value (as it happens instead for the subsequent 
fluid- mode oscillations). Section IV D is devoted to a 
thorough discussion of these issues; for the moment, we 
simply ignore this problem and focus our attention only 
on the part of the waveform dominated by fluid modes. 

Each panel of Fig. 6 compares the ID, exact 'I^'o 
(dashed lines) with that computed via the 3D code (solid 
lines) for the four values of the perturbation A. The ex- 
traction radius is (in both codes) f = 110 and this im- 
plies that a nonzero, constant offset for u < is present. 
Note, in this respect, the good consistency between 3D 
and ID results for u < 0, confirming here the informa- 
tion enclosed in the bottom panel of Fig. 2. After the 
initial high-frequency (unphysical) oscillations, the top- 
left panel of Fig. 6 shows that an excellent agreement 
between the waveforms is found when the perturbation 
is small. Then, for larger values of A (until it assumes val- 
ues that cannot be considered a perturbation anymore) 
the amplitude of the oscillation in the 3D simulations 
becomes smaller with respect to the linear case, suggest- 
ing that nonlinear couplings (specifically, couplings with 
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FIG. 6: (color online) Metric waveforms extracted at f = 110 computed from 3D simulations (dashed lines) and ID linear 
simulations (solid lines) for different values of the perturbation A. If A is sufficiently small {e.g. , A < 0.01) the outcomes of the 
two codes show good agreement. If the "perturbation" is large (A ~ 0.1), nonlinear effects become dominating. 



overtones as well as couplings with the radial modes) are 
redistributing the energy of the I = 2, m = oscillations 
triggered by the initial perturbation. In Sec. IV F, we will 
argue that couplings between modes become more and 
more relevant when the perturbation increases, giving a 
quantitative explanation to the phenomenology that we 
observe. This effect is summarized in Fig. 7, which dis- 
plays the amplitude A20 obtained by fitting the wave- 
form with the template Eq. (70) versus the magnitude 
of the perturbation for ID (linear) and 3D (nonlinear) 
simulations. It is evident from the figure that there is a 
consistent deviation from linearity already when the per- 
turbation is relatively small (A < 0.02). As a measure 
of the global agreement between ID and 3D waveforms 
(as a function of the initial perturbation A) we list in Ta- 
ble VI the I 2 residuals 1Z = 1 - 9(* 1D , * 3D ) and the l°° 
distances X>(^ 1D ,* 3D ). 



TABLE VI: "Global-agreement" measures computed on the 
interval Au — [50, 3000] (after the junk burst) at the outer- 
most detector. Here, TZ = 1 - 9(* 1D , * 3D ) is the I 2 residual 
while £>(* 1D , * :iD ) is the l°° distance. 



A 


n 


V 


AO 


3.07 x 10~ 2 


7.93 x 10" 5 


Al 


4.88 x 10~ 2 


7.79 x 10~ 4 


A 2 


1.63 x 10 _1 


2.78 x 10~ :i 


A3 


9.96 x 10" 1 


2.04 x 10~ 2 



The 3D waveforms for AO and Al turn out to be 
damped on a time scale of about 20 ms. This damping 
time is much shorter than the one of the / mode or pi 
mode, as computed via the ID approach. This effective- 
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TABLE VII: Frequency analysis of the 3D waveforms (see 
Fig. 6) over the interval u £ [50, 3000]. The frequencies from 
ID simulations are uj D = 1581 Hz and = 3724 Hz. From 
left to right the columns report the amplitude of the pertur- 
bation, the /-mode frequency, its relative difference with the 
ID value, the pi-mode frequency and its relative difference 
with the ID value. 



A 


vf 3 [Hz] 


Diff.[%] 


C [Hz] 


Diff.[%] 


AO 


1578 


0.2 


3705 


0.5 


Al 


1576 


0.3 


3705 


0.5 


A2 


1573 


0.5 


3635 


2.4 


A3 


1623 


2.7 


3565 


4.3 



viscosity damping time t visc (that is related to the inverse 
of the viscosity coefficient) can be extracted by means 
of the fit analysis discussed above for the waveform. We 
have found that r VISC depends on the initial perturbation, 
being r visc ~ 0.022,0.132,0.203,0.129 s respectively for 
A = AO, Al, A2, A3. The best agreement with the expected 
physical value of T20 = 0.268 s is obtained for A = A2; 
both for larger and smaller perturbations the 3D results 
show even shorter damping times. The errors on these 
quantities are of the order of 0.5%. The interpretation 
of these results may include two different effects. The 
smaller damping time of the wave for the A = A3 pertur- 
bation with respect to the A = A2 one may be interpreted 
as due to the nonlinear couplings that allow the disexci- 
tation of the fundamental mode in other channels; as it 
can be seen from Fig. 7 and Fig. 19, the importance of 
nonlinear effects is larger for the simulation with pertur- 
bation A = A3. However, for perturbations smaller than 
A = A2 the effective viscosity is not found to decrease to- 
wards the expected perturbative value, as it could have 
been expected from the above argument. This discrep- 
ancy might be due to the numerical viscosity proper of 
the evolution scheme. Such numerical viscosity would 
have a bigger influence in low-perturbation simulations, 
where the energy lost from the fundamental mode into 
other modes is smaller (while in higher-perturbation sim- 
ulations the coupling of modes is the dominant effect). 
Although the detailed analysis of the numerical viscosity 
of the 3D code is beyond the scope of the present work, 
we checked that, as expected, it depends on the grid reso- 
lution. We performed tests using a three-refinement-level 
setup with the resolution of the coarsest grid (with lim- 
its [—120, 120] in the three directions) set at the values 
^xyz — 2 (low), 1 (medium) and 0.5 (high). Using these 
three resolutions, we observed that, in the case of the 
coarsest grid, there was an initial explosion in the am- 
plitude, then followed by a strong damping during the 
first five gravitational-wave cycles. This shows that this 
resolution is not even sufficient to extract the qualita- 
tive behavior of the waveform. On the other hand, the 
other two resolutions did not show any qualitative dif- 
ference in addition to the different value of the "effective 
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FIG. 7: (color online) Comparison between the oscillation am- 
plitude A20 (from fits) in the ID (linear) and 3D (nonlinear) 
simulations versus the initial perturbation A. Deviations from 
linearity are occurring already for very small values of A. 



viscosity" , that is smaller for higher resolutions. We also 
checked whether there is a measurable effect due to the 
artificial atmosphere. Focusing only on the A = AO per- 
turbation, we varied the value of the rest-mass density 
of the atmosphere in the range p atm = 10[~ 5 ' _6 ' _7 ^p max , 
without finding any significant influence on the values of 
r vlsc . We leave to forthcoming studies a detailed analysis 
of the viscosity of the 3D evolution code. 

Finally, we have also Fourier-transformed the 3D wave- 
forms to extract the fluid-mode frequencies and we have 
compared them with the linear ones. This comparison is 
shown in Table VII. Apparently, the frequency of the / 
mode (that dominates the signal) is less sensitive to non- 
linear effects than its amplitude, as it can be seen from 
the fact that only the A = A3 initial data are such to force 
the star to oscillate at a frequency slightly different from 
that of the linear approximation. On the other hand, the 
first overtone (the pi mode) seems more sensitive. It is 
in any case remarkable that for A = AO and A = Al the 
frequencies from 3D and ID simulations coincide at bet- 
ter than 1%, suggesting that the main gravitational- wave 
frequencies are only mildly affected by nonlinearities. 



C. Nonradial oscillations: comparing ID and 3D 
curvature waveforms 

This section is devoted to the comparison between ID 
and 3D curvature waveforms. In the ID code one can use 
the relation 

rli m = rh^ = N t (¥& + i¥°l) (73) 

to obtain the Newman-Penrose scalar (multiplied by the 
extraction radius) rip^ m from the gauge-invariant metric 
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FIG. 8: (color online) ID versus 3D evolution of the rip%° 
curvature waveform for perturbation A = AO. The initial 
transient is consistent between the two evolutions. The in- 
set concentrates on the initial part of the waveform. 
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FIG. 9: (color online) Comparison of the PSD of various wave- 
forms. The PSD of the ripf and from the ID code are 
compared with those from the 3D code, obtained after inte- 
gration. The initial perturbation is A = AO. The in-mode 
frequencies W\ and W2 are superposed to the spectra as ver- 
tical dot-dashed lines. See text for discussions. 



master functions. Because of our choice of initial con- 
ditions, we shall consider only "F^o m the following. 5 
The second time-derivative of is computed via fi- 
nite differencing, by applying twice a first-order deriva- 
tive operator with 4th-order accuracy. By contrast, in 
the 3D code ipf 3 is extracted independently of the met- 
ric waveform. Then, one computes rtpf 1 , where r is an 
approximated radius 6 from Eq. (59) with M = 1.4. 

Figure 8 displays the r^|° waveforms from ID (solid 
line) and 3D (dashed line) evolutions with perturbation 
A = AO. The extraction radius is f = 110 in both codes. 
Visual inspection of the figure immediately suggests that: 
(i) The initial transient in the ID metric waveform pre- 
ceding the setting in of the quasiharmonic /-mode oscil- 
lation results in a highly damped, high-frequency oscil- 
lation; (ii) The initial transient radiation has the same 
qualitative shape in both the ID and 3D waveforms, al- 
though the amplitude of the oscillation is larger in the 
latter case. At this point one should note that: (i) In 
the ID case, although the conformally flat condition is 
imposed at t = 0, the constraint is solved numerically 
and thus a small violation of this condition occurs; (ii) 
The violation is expected to be larger in the 3D case, be- 
cause of the larger truncation errors. It is in any case re- 



5 Note that in principle one could compute 1/14 independently, solv- 
ing the Bardeen-Press-Teukolsky equation [92]. 

6 This is an approximate relation as f is a coordinate radius and the 
mass inside the sphere of radius f is time dependent. We neglect 
all higher-order effects here as this approximation is sufficiently 
accurate for our purposes. 



markable that, as the figure shows, these errors (e.g. the 
slightly different shapes of /C20, the linear interpolation 
from spherical to Cartesian coordinates, etc.) are suf- 
ficiently under control to produce the same qualitative 
behavior besides small quantitative differences in the ini- 
tial part of the ID and 3D waveforms. 

The question that occurs at this point is whether the 
violation of the conformally flat condition introduces 
some amount of physical w-mode excitation in the wave- 
forms. To answer this question we show in Fig. 9 the 
Fourier power spectral density (PSD) of the rijff wave- 
forms of Fig. 8. The PSD is computed all over the wave- 
form and not only during the "ring-down" , because of the 
difficulty of separating reliably this part from the "pre- 
cursor" [10]. We are aware of the problems related to 
the precise determination the u>-modc frequencies and to 
their location in the waveform (see Refs. [10, 93] for a 
related discussion), and in particular of the fact that the 
Fourier analysis can not provide accurate and definitive 
answers, essentially because, in the presence of damped 
signals with frequencies comparable to the inverse of the 
damping time, the Fourier spectrum results in a broad 
peak. However it represents a fundamental part of the 
analysis and, in the present case, is preferable to a fitting 
procedure because of the already mentioned problem of 
separating the precursor from the ring-down part. 

The dashed-dotted vertical lines of Fig. 9 locate the 
first two w-mode frequencies of this model, v Wl = 10.09 
kHz and v W2 = 17.84 kHz. These frequencies have been 
computed by K. Kokkotas and N. Stergioulas via an in- 
dependent frequency-domain code and have been kindly 
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given to us for this specific comparison. Two of the max- 
ima of the PSD of the 3D rip 20 waveform can be asso- 
ciated to the frequencies v Wl and v W2 , even if they are 
in a region very close to the noise. The frequency u Wl is 
probably slightly excited also in the ID case (see inset), 
while only noise is present around v W2 . In Fig. 9 we show 
also the PSD of the ID ty^o (dashed dark line) and the 
one of the 3D $20 (dashed light line) obtained from the 
double time-integration of ip 20 . In both cases, it is not 
possible to disentangle v Wx and is W2 from the background 
noise. 

The fact that a signal characterized by highly damped 
modes is much less evident in the PSD of the metric wave- 
form than in the corresponding curvature one is simply 
due to the second derivative that relates the two gauge- 
invariant functions. When space-time modes are excited, 
the metric waveform is (approximately) composed by a 
pure ring-down part plus a tail contribution [94] , that is 

^20 * e ~ at + Pt~ 7 > where a = a + iuo (a is the inverse 
of the damping time and u> the w-mode frequency) and 
(3 is a numerical coefficient. When one takes two time 
derivatives to compute rip 20 from vE^c/ > the tail contribu- 
tion is suppressed by a factor t~ 2 and the oscillatory part 
of the waveform emerges more sharply. This comparison 
suggests that the best way to extract information about 
w modes (especially when their contribution is small) is, 
in general, to look at ripl m . In addition, it also highlights 
that, while it is not possible to exclude the presence of 
w modes in the ripf™ signal due to the small violation 
of the conformally flat condition at t = 0, at the same 
time we can not definitely demonstrate that those high 
frequencies present near the noise are attributable to w 
modes. In the next section we are going to show similar 
analyses on the spectra computed from Vf^o waveforms 
extracted a la Abrahams-Price from the 3D simulation. 

Finally, the global-agreement measures on the ipi ex- 
traction are K ~ 1.42 x 1CT 2 and V ~ 9.09 x 1(T 7 , and 
they highlight some differences between the linear and 
the nonlinear approach. 

The analysis discussed so far indicates that, in the 
present framework, the wave-extraction procedure based 
on the Ncwman-Pcnrose scalar ^4 seems to produce 
waveforms that, especially at early times, are more ac- 
curate than the corresponding ones extracted via the 
Abrahams-Price metric-perturbation approach. How- 
ever, one of the big advantages of the latter method 
is that the waveforms h + and h x are directly available 
at the end of the computation, and thus ready to be 
injected in some gravitational-wave-data-analysis proce- 
dure. By contrast, if we prefer to use Newman-Penrose 
wave-extraction procedures (which are the most common 
tools employed in numerical-relativity simulations nowa- 
days), we must consistently give prescriptions to obtain 

\B^ o) fromVl™- To do so, one needs to perform a double 
(numerical) time integration, with at least two free inte- 
gration constants to be determined to correctly represent 
the physics of the system. Inverting Eq. (73) following 




FIG. 10: (color online) Recovery of M/go fr° m two successive 
time integrations of ripf 3 extracted from the 3D simulation 
with initial perturbation A = AO. After the subtraction of a 
quadratic floor the waveform correctly oscillates around zero. 

the considerations of Sec. II C 1, we obtain the following 
result [see Eq. (38)]: 

rh«™ = N e (¥£ l + i¥£) (74) 

= / dt' [ dt"ripj m (t") + Q + Q 1 t + Q 2 t 2 , 
Jo Jo 

= rh im (t) + Q + Q 1 t + Q 2 t 2 , 

where Qo, Q\ and Q 2 are (still) undetermined integra- 
tion constants, which are complex if m 7^ 0. Note that 
this relation does not involve the Q 2 integration constant 
only if finite-radius extraction effects can be considered 
negligible (see below). 

Our aim is to recover the metric waveform that cor- 
responds to the 3D rip 20 waveform that we have charac- 
terized above. We consider the waveform of Fig. 8 up to 
t = 1500, where the reduction in amplitude due to numer- 
ical viscosity is already of the order of 30% with respect 
to the exact linear waveform. This sampled curvature 
waveform is integrated twice in time, from t = without 
fixing any integration constant to obtain rh£ m (t). 

The raw result of this double integration is shown in 
Fig. 10. The "average" of the oscillation does not lay 
on a straight line, as it does instead in the case of the 
waveforms of binary black-hole coalescence discussed in 
Ref. [64], but rather it shows also a quadratic correc- 
tion due to the finite extraction radius (see discussion in 
Sec. II CI). 

Indeed, when a "floor" of the form P(t) = Qo + 
Qit+Q 2 t 2 is subtracted, the resulting metric waveform is 
found to oscillate around zero, as it can be seen in Fig. 10 
and Fig. 11, which focuses on the beginning of the oscilla- 
tion. The values of the coefficients of P(i) obtained from 
the fit are Q = -4.338 x 10~ 7 , Qi = -1.2462 x 10" 7 
and Q 2 = -6.2046 x 10~ 9 . The fact that Q < is con- 
nected to the choice of initial data we made (i.e. k 2 Q 7^ 
at t = 0). Then, Qi 7^ indicates that the system is 
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FIG. 11: (color online) Selecting the best fitting function 
for the 3D integrated-and-subtracted metric waveform. The 
dash-dotted line refers to the curve obtained with the subtrac- 
tion of a cubic fit and the dashed line to the curve obtained 
with the subtraction of a quadratic fit. They are very similar 
to each other and to the ID waveform (solid line). Contrary 
to the case of the extracted metric waveform, no initial burst 
of radiation is present in the 3D integrated-and-subtracted 
waveform. 



(slightly) out of equilibrium already at t = and it is 

thus emitting gravitational waves since (0) ^= 0. This 
is consistent with the choice of initial data we made, that 
is a perturbation that appears instantaneously at t = 
without any radiative field obtained from the solution of 
the momentum constraint (since we use time-symmetric 
initial perturbations, for which fc20 = X20 = 0). 

We tested the robustness of the quadratic fit by adding 
a cubic term Q^fi to P(t) and then fitting again. In 

Fig. 11, we compare the 3D waveform corrected 

with a cubic fit (dashed line) with the one corrected with 
a quadratic fit (dash-dotted line) and with the "exact" 
ID metric waveform (solid line) output by the PerBACCo 
code. Note that the ID waveform has been suitably 
timcshiftcd in order to be visually in phase with the 
others at the beginning of the simulation. The figure 
suggests that the effect of the cubic correction is almost 
negligible (one only finds slight changes in the very early 
part of the waveform). The values of the fitting coeffi- 
cients Q t are Q = -1.0096x10 A Qi = -6.3331x10 A 
Q 2 = -4.747 x 10" 8 , Q 3 = 6.7114 x 10~ 14 . The fact that 
Q3 is many orders of magnitude smaller than the other 
coefficients is a good indication that the quadratic be- 
havior is indeed the best choice here. Consistently with 
the curvature waveform of Fig. 8, we note the excellent 
agreement between ID and 3D (integrated) metric wave- 
forms also in the initial part of the waveform, i.e. up 
to t ~ 200 (corresponding to the high-frequency oscil- 
lation in rip1°). Evidently, this is in contrast with the 
Abrahams-Price metric waveform in the top-left panel of 
Fig. 6 (we will elaborate more on this in the next section). 



Finally, we point out that the coefficient Q2(r) shows, 
as expected, a clear trend towards zero for increasing 
values of the extraction radius. 



D. Advantages and disadvantages of the 
Abrahams-Price metric wave-extraction procedure 

The analysis carried out so far suggests that both the 
Regge- Wheeler- Zerilli metric-based and the Newman- 
Penrose ^4-curvature-based wave-extraction techniques 
can be employed to extract reliable gravitational wave- 
forms from simulations of compact self-gravitating sys- 
tems. For the particular case of an oscillating neutron 
star as considered here, both extraction methods allow 
to obtain waveforms that are in very good agreement 
with the linear results. Despite this success, the two ap- 
proaches are not free from drawbacks. Let us first focus 
on r"0| m curvature waveforms. The comparison between 
ID and 3D r^f waveforms in Fig. 8 (as well as between 
integrated metric waveforms in Fig. 11) shows good con- 
sistency between the two (as long as the effects of nu- 
merical viscosity on the evolution of the system remain 
negligible). As we mentioned above, we think that the 
most important information enclosed in Fig. 8 is that 
the differences between the high-frequency oscillations in 
the initial part of the waveforms (where w modes are 
probably present in the 3D case) are small. This fact 
makes us confident that the violation of the 3D Hamil- 
tonian constraint at t = (due to its approximate so- 
lution 7 ) as well as the violation of the conformally flat 
condition are sufficiently small to avoid pathological be- 
havior during evolution. A further confirmation of the 
accuracy of the evolution and of the curvature extrac- 
tion is given by Fig. 12: The quantities ripf 1 extracted 
at various radii (f € {30,60,90,110}) and plotted ver- 
sus retarded time are all superposed. This confirms the 
theoretical expectations of the peeling theorem [95] and 
indicates (once more) that the quantity r0|° is accu- 
rately computed. In Fig. 12, r is obtained from r via 
Eq. (59). The retarded time is approximated with the 
standard r*, where the constant mass M = 1.4 has been 
used. In our setup, the only subtle issue about ripi™ 1 
seems to be the computation of the corresponding met- 
ric waveform via a double time integration. Although 
we were able to obtain a rather accurate metric wave- 
form, the time-integration procedure (including the eval- 
uation of the integration constants) may not be likewise 
straightforward in other physical settings. By contrast, 
the Abrahams-Price wave-extraction procedure directly 
produces the metric waveform and no time integrations 



7 We recall that the 3D Hamiltonian constraint is solved at the 
linearized level on an isotropic grid and then the resulting met- 
ric perturbation is interpolated on the Cartesian grid. Typi- 
cally, this procedure leads to larger errors than if solving the 
constraints directly on the Cartesian grid. 
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FIG. 12: (color online) The quantities ripl extracted at dif- 
ferent radii are superposable (as expected from the "peeling" 
theorem). See text for further explanation. 



are needed. For this reason, it looks a priori more ap- 
pealing than if; 4 extraction. Unfortunately, the results 
that we have presented so far (notably our Fig. 6) indi- 
cate that this computation can be very delicate and can 
give unphysical results even in a very simple system like 
an oscillating polytropic star: we have found that 
extracted in this way is unreliable at early times, because 
of the presence of high-frequency, highly damped oscil- 
lations, that are instead absent in both the ID linear 
metric waveforms and the 3D metric waveforms time- 
integrated from n/>f . The unphysicalncss of this initial 
"burst" of radiation is evident from Fig. 13, where the 
extractions at various radii f € {30,60,90,110} of the 

quantity ^^o' are compared: the amplitude grows with 
f, instead of decreasing progressively to approach a con- 
stant value (as it is the case for the /-mode-dominated 
subsequent part of the waveform) 8 . The weird behavior 

at early times of the extracted indicates that this 
function does not satisfy the Zerilli equation in vacuum. 
Consistently, the perturbative Hamiltonian constraint in 
vacuum, Eq. (64) with H 2 o = 0, constructed from the 
3D metric multipoles (x20j ^20 )j must be violated of some 
amount in correspondence of the junk 9 . This reasoning 
suggests that the junk may be the macroscopic manifes- 
tation of the inaccuracy in the initial-data setup at t = 
(i.e. of solving the linearized Hamiltonian constraint first 



8 To assess this statement we have also performed simulations with 
extraction radii up to f = 200. 

9 The Zerilli equation, and thus the Zcrilli-Moncricf master func- 
tion, is obtained by combining together the perturbative Einstein 
equations, one of which is precisely the perturbative Hamiltonian 
constraint in vacuum. The Zerilli equation is satisfied if and only 
if the perturbative Hamiltonian constraint is satisfied too. See 
for example Ref. [16] for details. 
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FIG. 13: Metric waveforms computed via the Abrahams-Price 
procedure and extracted at different radii, for a 3D evolution 
with perturbation A = AO. Note how the amplitude of the 
initial burst grows linearly with f. 



and then interpolating) , possibly further amplified by the 
wave-extraction procedure. This statement in itself looks 
confusing, because we have learned, from the analysis of 
■04, that the Einstein (and matter) equations are accu- 
rately solved and that the errors made around t = due 
to the violation of Hamiltonian constraint are relatively 
negligible. The relevant question is then: is it possible 
that small numerical errors, almost negligible in rip 20 , 
may be amplified in at such a big level to produce 
totally nonsensical results? The following discussion pro- 
poses some heuristic explanation. 

To clarify the setup of our reasoning, let us first remind 
the reader of the basic elements of the Abrahams-Price 
metric wave-extraction procedure and, in particular, the 
role of Eq. (4). At a certain evolution time t, the nu- 
merical metric g^vit) is known at a certain finite accu- 
racy on the Cartesian grid. One selects coordinate ex- 
traction spheres of coordinate radius f = \/ x 2 + y 2 + z 2 
on which the metric is interpolated via a second-order 
Lagrangian interpolation. Isotropic-coordinate systems 
(f, 9, 4>) naturally live on these spheres and thus one de- 
fines spherical harmonics. Then, the metric <7 M „ is for- 
mally decomposed in a Schwarzschild "background" g ^ 
plus a perturbation h^. The next step is to choose a 
coordinate system in which the background metric is ex- 
pressed. The standard approach is to use Schwarzschild 
coordinates, although this choice actually introduces sys- 
tematic errors that may relevantly affect the waveforms. 
This has been recently demonstrated in Ref. [6]. Al- 
though we are aware of this fact, we prefer to neglect 
this source of error, on which we will further comment 
below. Choosing Schwarzschild coordinates means that 
one needs to compute a Schwarzschild radius r. This is 
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given by the areal radius of the extraction two-spheres. 
Proceeding further, h^ v is decomposed into seven (gauge- 
dependent) even-parity (f/o, Hi, H 2 , , G, K) and 
three (gauge-dependent) odd-parity multipoles (that we 
do not consider here). From combinations of the seven 
even-parity multipoles and of their radial derivatives, see 
Eqs. (41) and (42) of Ref. [19], one obtains the gauge- 
invariant functions ki m and Xlm% as well as the derivative 
d r ke m . The last step is the computation of the Zerilli- 
Moncrief function via Eq. (4). Various sources of errors 
are present. In particular, we mention the errors origi- 
nating from: (i) the discretization of g^ (and its deriva- 
tives), from the numerical solution of Einstein's equa- 
tions; (ii) the interpolation from the Cartesian grid to the 
isotropic grid; (hi) the computation of the metric multi- 
poles via numerical integration over coordinate (gauge- 
dependent) two-spheres. Our aim is to investigate how 
these inaccuracies on (xim,kt m ,d r ki m ) can show up in 
^hn a t l ar § e extraction radii. In the limit r 3> M, Eq. (4) 
reads 
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where Zi m = Xim ~ rd r ki m + A/2ki m and r is the areal 
radius of the coordinate two-spheres. The Abrahams- 
Price wave-extraction procedure introduces then errors 
both on r and Zg m . In particular, the errors on the 
(gauge-invariant) multipoles (xim,kg m ,d r kg m ) conspire 
in a global error on Zg m . In a numerical simulation one 
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tions of the perturbation equation on a Schwarzschild 
background, and r w is the radial Schwarzschild coor- 
dinates; 5Zi m encompasses all possible errors due to the 
multipolar decomposition procedure, and Sr various inac- 
curacies related to the determination of the areal radius 
(e.g. , those related to gauge effects). As a result, for the 
"extracted" Zerilli-Moncricf function we can write 
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This equation shows that, if SZi m is not zero at a certain 
time (and does not decrease in time like l/ r Schw ) there is 
a contribution to the global error on that grows lin- 
early with the extraction radius. This qualitative picture 
is consistent with what we observe in the 3D waveforms: 
a small error on SZ20 introduced at t = 0, because of the 
approximate solution of the constraints (as indicated by 
the analysis of rip1° curvature waveforms), can show up 
as a burst of radiation whose amplitude increases linearly 
with the observer location. Note that what really counts 
here is the error budget at the level of (X20, k 2 $, d r k 20 ) 
and the related violation of the pcrturbativc Hamiltonian 
constraint, Eq. (64). Indeed, it might occur that, even 



if the 3-metric 7^ is very accurate and the constraints 
are well satisfied at this level, the extraction procedures 
adds other errors (for example due to the multipolar de- 
composition, computation of derivatives etc.) that may 
be eventually dominating in SZ 2 q ■ This observation may 
partially justify why r"0|° is well behaved, while ty^o 1S 
not. Finally, we note that in our evolution 5r is typi- 
cally very small, so that we have r Schw r* r with good 
accuracy. 

Because of the complexity of the 3D wave-extraction 
algorithm, we were able neither to push forward our level 
of understanding, nor to precisely diagnose the cause of 
the aforementioned errors 10 . This is now beyond the 
scope of the present work and will deserve more atten- 
tion in the future. By contrast, we can exploit the simpler 
computational framework offered by the ID PerBACCo 
code to "tune" the error SZi rn in order to produce some 
initial "spurious" burst of radiation, and then possibly 
observe that its amplitude grows linearly with f. In 
the ID code Sr is zero by construction, so that all er- 
rors are concentrated on SZi rn . The constrained scheme 
adopted in the pcrturbativc code (which is second-order 
convergent) allows to accurately compute the multipoles 
(Xim, kim) at every time step, and the Hamiltonian con- 
straint is satisfied by construction. Then, d r ke m is ob- 
tained via direct numerical differentiation of kg m . Con- 
sequently, the error SZi m depends on the resolution Ar 
as well as on the order of the finite-differencing represen- 
tation of d r ki m . 

In the following we shall analyze separately the effect of 
resolution and of the approximation scheme adopted for 
the numerical derivatives. First, we approximate d r kg m 
with its standard first-order finitc-diffcrcncing represen- 
tation, i.e. d r ki m f» (fcjlp! — kj)/Ar and we study the 

behavior of the extracted computed using Eq. (4), 

versus extraction radius and resolution. Second, we use 
a fixed Ar, but we vary the accuracy of the finite- 
differencing representation of d r ke m , contrasting first- 
order, second-order and fourth-order stencils. The results 
of these two analyses, for I = 2, m — 0, are shown in 
Figs. 14 and 15 respectively. In the top panel of Fig. 14 
d r k 2 Q is approximated at first-order, with a resolution 
of 10 points inside the star (Ar ~ 0.9). This resolu- 
tion approximately corresponds to the resolution of the 
coarsest refinement level used in the 3D code. The ex- 
tracted waveform ty^d is shown at different observers, 
f G {30,60,90,110}: An initial burst of junk radia- 
tion develops at early times and its amplitude grows lin- 
early with the extraction radius (and it keeps growing for 
f > 110). This behavior looks identical to that found in 
the 3D simulations. In the bottom panel of Fig. 14 we fo- 
cus on f = 110 only, but vary the number of radial points 



10 For example, we mention, in passing, that we have also tried 4th- 
order Lagrangian interpolation, without any visible improvement 
on the waveform. 
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FIG. 14: (color online) Metric waveforms from linear ID evo- 
lutions. Top panel: Low resolution simulation. A burst 
of junk radiation at early times is present and its amplitude 
grows linearly with the extraction radius. Bottom panel: 
By increasing the resolution, the initial junk disappears. 



N r inside the star, namely N r e {10,20,30,40,50, 100}. 
The figure shows that the initial junk is not present at 
higher resolutions (N r > 20) and that the waveform con- 
verges to the exact profile 11 . But varying the resolution 
only shifts the occurrence of the burst at farther radii: 
observers at f 3> 110 still see this burst appear and grow 
linearly with f . 
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11 As discussed in Rcf. [10] we cross-checked things also by matching 
the Zcrilli-Moncrief function at the surface and evolving it with 
the Zerilli equation outwards. We found good agreement between 
the "matched" and the "computed" waveforms. 
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FIG. 15: The amount of junk radiation present in the wave- 
forms extracted in the ID code is smaller when higher-order 
differential operators are used in Eq. (4) to compute the the 
Zerilli-Moncrief function. 



A complementary analysis is shown in Fig. 15, where 
we fix the resolution at N r — 10 (for f = 110), but we 
change the accuracy of the numerical derivative d r k2o- 
As expected, the initial junk disappears when the accu- 
racy of the numerical differential operator is increased: a 
second-order operator produces only a small amplitude 
bump, that is not present when the fourth-order operator 
is employed. At this stage, the conclusion is clear: the 
convergence of the Zerilli-Moncrief function computed 
from the separate knowledge of the multipolcs kg m and 
Xtm is a delicate issue that must be analyzed with care 
according to the physical problem under consideration. 
The violation of the perturbative Hamiltonian constraint 
and, in particular, the accuracy of the numerical deriva- 
tive d r ki m (note that we refer to the induced violation 
at the level of the wave extraction and not at that of 
the solution of the perturbation equations) seems to play 
an important role in the convergence properties of the 
waveforms. The main conclusions of the aforementioned 
numerical tests are: (i) The errors in seem to behave 
like suggested in Eq. (77); (ii) The phenomenon occurs 
in the same way in both the ID and 3D code, although 
the fine details of the oscillation arc different. 

Focusing on the ID PerBACCo code, an accurate 
is obtained using sufficiently high resolution (N r = 300) 
as well as a fourth-order representation for d r ki m . These 
prescriptions are accurate enough for the problem ad- 
dressed in this work, although they may not be sufficient 
for other stellar models or other initial perturbations. For 
example, using the PerBACCo code, with the same initial 
data setup discussed here, in order to study the time evo- 
lution of perturbations of stars with realistic EOS proved 
that higher resolutions are typically needed to produce 
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FIG. 16: (color online) PSD of the 3D metric waveforms 
(dashed line) extracted a la Abrahams-Price (at f = 110) 
and the corresponding ID waveform (solid line) for a simula- 
tion with perturbation A = AO. The Fourier spectrum of the 
junk radiation is compatible with some lii-mode frequencies. 
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FIG. 17: (color online) The top panel compares the standard 
Zerilli-Moncrief $20 (depicted with lines) and the generalized 
^ r |J MP (depicted with point markers) for the first part of the 
gravitational-wave signal in a simulation with perturbation 
A = Al. The waves extracted at four radii f £ {30, 60, 90, 110} 
are shown. The bottom panel shows that the present junk 
radiation negligibly depends on the choice of the coordinates 
of the background metric. 



convergent waveforms of comparable accuracy [11]. Like- 
wise, for a polytropic EOS and initial data given by a 
Gaussian pulse in 'J/^'o > the same Ref. [11] showed that at 
least fourth-order accuracy in d r k2o is needed in order to 
have a consistent extraction of already at t = (see 
Appendix A of Rcf. [11]). This suggests that the presence 
of linearly growing junk radiation in the computation of 

^tm fr° m the multipoles {kt m ,xim) can appear ubiqui- 
tously in the time evolutions of the perturbation equa- 
tions with the PerBACCo code. The presence of this junk 
m 

*/m is ttle 

macroscopic manifestation of the violation 
of the pcrturbative Hamiltonian constraint due to errors 
(notably, in the discretization of the derivatives) intro- 
duced in the wave-extraction procedure. These (typically 
small) numerical errors are eventually magnified by the 
presence of an overall r factor in Eq. (75). Note that this 
phenomenon occurs even if the computation of the mul- 
tipoles {ki m ,xim) is very accurate and the Hamiltonian 
constraint is satisfied by construction in the evolution 
algorithm. The analysis that we have presented here 
suggests that either increasing the resolution or, more 
reasonably, implementing higher-order differential oper- 
ators in the pcrturbative "extraction" procedure are vi- 
able proposals to compute convergent waveforms. 

In the 3D case the situation is more involved and we 
have not succeeded in making statements as solid as in 
the ID case. We can only rely on analogies: (i) The 
appearance of the junk occurs in a way similar to the 
ID case when the accuracy of the ID Zerilli function is 
low; (ii) The two time evolutions look qualitatively very 
similar. Yet, it is not technically possible to use in the 



3D code resolutions equivalent to those of the ID code. 
By analogy with our perturbative results, we can only 
conclude that it is not unreasonable that the junk in the 
3D waveforms is the macroscopic manifestation of inac- 
curacies hidden in the implementation of the Abrahams- 
Price wave-extraction procedure. The analysis presented 
here points out that such metric wave-extraction proce- 
dures require typically more subtle care than expected 
and these subtleties must be kept in mind in developing 
more modern wave-extraction routines. 

The PSD of in both the ID case (solid line) and 3D 
case (dashed line) is displayed in Fig. 16. The perturba- 
tion is A = AO and the extraction radius is f = 110. The 
spectrum of the 3D Zerilli waveform is consistent with 
what we observed in ripf below 10 kHz, but it looks dif- 
ferent at higher frequencies (compare it with that of the 

"integrated" 'F^'o m Fig- 9): here the PSD shows broad 
peaks attributable to the initial part of the waveform. 
Recovering the reasoning started in the previous section 
about the presence of the w modes, we can observe that 
the frequencies contained in the junk are compatible with 
the w-modc frequencies u Wl and v W2 (indicated by dashed 
vertical lines in Fig. 16). However, since such frequencies 
belong to an unphysical part of the waveform, we prefer 
to consider them unphysical as well. 

We conclude by mentioning, in passing, that the ini- 
tial junk radiation is essentially not related to the sys- 
tematic error introduced by fixing Schwarzschild coor- 
dinates for the background metric g® v . This fact is sug- 
gested by Fig. 17, where we contrast the standard Zerilli- 
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Moncrief function (which assumes Schwarzschild 

coordinates for the background) with the generalized 
*2™ P one based on the Sarbach-Tiglio [55] and Martel- 
Poisson [56] perturbation formalism, which does not re- 
quire any gauge-fixing condition for the background sub- 
manifold A'/ 2 . This particular simulation was performed 
over a grid with three refinement levels and cubic boxes 
with limits [-120, 120], [-24,24] and [-12,12]. The resolu- 
tion of each box is coarser than in the previous simula- 
tions, namely A xyz = 1.875, 0.9375 and 0.46875 respec- 
tively Evidently with this resolution the waveforms are 
less accurate, but we do not mind at this stage, since we 
are interested in an intrinsic comparison between extrac- 
tion procedures at fixed resolution. The function ^fZ MP 
that we use is given by the straightforward 12 implemen- 
tation of Eq. (4.23) of Ref. [56]. Note that this expression 
is equivalent to the combination of Eqs. (20), (25), (26) 
and (27) of Ref. [55]. The top panel of Fig 17 displays 

^20 (lines) and 'I'fo" 1 ^ 15 (point markers) for observers at 
f g {30,60,90,110}. It highlights that the differences 
in the early-time part of the waveforms are very small. 
By contrast, the bottom panel of the figure, showing the 
difference — t J , 2™ P j indicates that removing (part 
of) the systematic errors generates some improvement, 
but this is too small to be of any relevance. This analy- 
sis suggests that the inaccuracies in the early-time part 
of the waveform are essentially not related to the spe- 
cific computation of the (generalized) Zerilli function, 
but rather connected to the underlying multipolar ex- 
traction infrastructure (grid setup, approximate solution 
of the constraints, interpolation procedures, computation 
of the derivatives of the metric etc.), on which we have 
relatively little control. A comprehensive analysis of the 
problems related to the generalized extraction procedure 
will be presented elsewhere [71]. We remark, however, 
that systematic effects that are very small in our physical 
system, as emphasized by our Fig. 17, may be not small 
in other situations, as found in Ref. [6]. For this reason, 
we emphasize that the formalism of Rcfs. [55, 56] is the 
actual correct metric formalism to extract waveforms out 
of a numerical space-time that can be considered a small 
deformation of the Schwarzschild one. As such, it must 
be taken into account properly in numerical codes. 



E. Generalized quadrupole-type formulas 

Finally, we study the performances of the various gen- 
eralized quadrupole-type formulas that we have intro- 
duced in Sec. II C. The results of our analysis are shown 



in Fig. 18. The left panel of the figure displays rh 20 wave- 
forms obtained via the SQF1, SQF2, SQF3 and SQF4 
[see Eqs. (52-55)] for perturbation A = AO. The right 
panel complements this information by showing (for sev- 
eral perturbation magnitudes A) the relative difference in 
amplitude between the various SQFs and the correspond- 
ing gauge-invariant Zerilli-Moncrief function \&20 ■ This 
analysis highlights that the quadrupolc formula gives an 
excellent approximation to the phasing of the actual sig- 
nals. By contrast, there is a systematic over- or under- 
estimation of the amplitude depending on the choice of 
SQF. 

A related observation is that the discrepancy between 
the quadrupolc formula and the gauge-invariant wave- 
form is not due to the fact that waveforms are extracted 
at a finite radius. Our results are consistent with those 
of Shibata and Sckiguchi [5] , who performed an analysis 
similar to ours (and also considered uniformly rotating 
stars) , but without the possibility of contrasting their re- 
sults with linear evolutions. In this respect, the main 
conclusion of Ref. [5] was that, although the amplitude 
of the waveform is systematically underestimated by the 
quadrupole formula (47), it is however sufficiently accu- 
rate to capture both the frequency and the phasing (that 
are the most important quantities for detection) of the 
waveforms in a proper way. These results are fully con- 
firmed here using totally different codes. 



F. Nonlinearities 

In this section we comment on the onset of nonlinear ef- 
fects for high values of the initial perturbation amplitude 
A, showing, for the first time using full GR simulations, 
evidences for mode couplings and for the appearance of 
nonlinear harmonics. 

In the linear regime (A = AO) the star is oscillating at, 
essentially, the frequency of the fundamental quadrupolar 
proper fluid (quasinormal) mode of pulsation. The prin- 
cipal linear modes excited are thus the (I, m) = (2, 0) 
one and its overtones. For growing values of the initial 
perturbation, we observe that in 3D simulations, differ- 
ently from the linear ones, the amplitude of the multipole 
does not increase proportionally to A but, instead, 
is progressively reduced (see Fig. 7). This fact could be 
interpreted as the results of a typical phenomenon in non- 
linear systems in which linear modes couple, generating 
nonlinear harmonics. Naively, one could think that the 
"energy" associated to the I = 2 mode is redistributed 
to the others while the system departs from the linear 
regime 13 . As we saw in Sec. IV A, radial modes of os- 



2 With this we mean that we do not take into account any time 
dependence of the background metric due to coordinate effects. 
This possibility can be anyway easily taken into account by the 
formalism. We postpone to a future work the related discus- 
sion [71]. 



We would like to stress that this picture is not so simple and 
rigorous, in particular there is no proof of the completeness of the 
star quasinormal modes (even in the nonrotating case) and the 
definition of an energy per mode is definitely not straightforward. 
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FIG. 18: (color online) Comparison between different wave-extraction procedures in the 3D simulations, in particular between 
the SQFs and the gauge-invariant wave-extraction procedure. The left panel shows the waveforms. The frequencies almost 
coincide, while the amplitudes are strongly underestimated in all cases except SQF1, where the amplitude is strongly overes- 
timated. The right panel shows the relative difference between the amplitudes of the various SQFs and the Abrahams-Price 
wave extraction. These amplitudes are estimated with the fit procedures. 
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FIG. 19: (color online) PSD of the quantity (p)2o{t) [see 
Eq. 78] for different values of the initial perturbation ampli- 
tude A. The spectra of the (I = 2, m = 0) mode obtained in 
simulations with larger perturbations contain more frequen- 
cies, which originate from the nonlinear couplings with the 
overtones and with the radial modes. 



cillation are already present in the evolution of the equi- 
librium models. As a consequence, we expect to reveal 
couplings between nonradial and radial modes (F and 
its overtones Hi,H2, ■■■) as a result of the onset of some 
nonlinear effect. In addition, we detect signals in multi- 



See the discussion in Ref. [93]. 



poles of 'J?^ with I = 4, 6 and m = 0, 4, that are even- 
parity axisymmetric modes and nonaxisymmetric modes 
triggered by the Cartesian grid. The amplitudes are very 
weak compared to those of the / mode for any value of A, 
typically 2 orders of magnitude smaller for £ = 4, m = 
and 3 orders of magnitude for the others, but in principle 
they are present and must be considered. As far as the 
odd-parity modes with m = 1, 2, 3, and t = 3, 5 are con- 
cerned, they are all forbidden by the symmetry imposed 
on the computational domain (octant). 

As a strategy to study nonlinearitics, we consider the 
rest-mass-density projections: 

(p)Ut) = J d 3 xp(t,x)Y; m (78) 

and we apply to them the Fourier analysis. Like all the 
global variables, p contains all the frequencies of the sys- 
tem. Its projections in Eq. (78) allow to separate the con- 
tribution of each mode (£,m). Fig. 19 shows the power 
spectrum of (p}20 for the four different values of A. The 
signal for A = AO contains the 3 frequencies of the lin- 
ear modes /, p\ and p2- The same happens for A = Al 
and A = A2: the amplitudes of the linear modes grow 
linearly with A and some new frequencies are present 
with small power for A = A2. In the case of A = A3, 
which corresponds to a pressure perturbation of 10% of 
the central TOV value, the spectra is rich of nonlinear 
harmonics. Most of them can be recognized as due to 
weak couplings, i.e. sums and differences of linear mode 
frequencies v\ ± i>2, also called combination tones. In 
particular we identify the nonlinear harmonics of the / 
mode and its overtone fipi and many frequencies f±F, 
f ± H i and p\ ± H\ due to the radial and nonradial mode 
couplings. Such couplings have been previously and ex- 
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tcnsivcly studied in Refs. [86, 89, 90] using the Cowling 
approximation as well as the conformally flat approxima- 
tion to GR. In addition, the couplings between radial and 
nonradial modes have been studied in detail in Ref. [96] 
by means of a second-order perturbative approach. Note 
how our fully general-relativistic results are consistent 
with all these studies. 

The projection (p)oo describes essentially the radial 
mode of pulsations; analyzing this quantity instead of 
{p)2o gives analogous results in term of couplings. From 
the analysis of higher multipoles we compute the frequen- 
cies of the linear modes, finding v = 2404 Hz for £ = 4 
and v = 2988 Hz for £ = 6. No couplings can be clearly 
recognized in these data. We stress that frequencies of 
the nonaxisymmetric modes (m = 4) are the same as the 
axisymmctric ones because the star is nonrotating and 
modes arc degenerate in m. 



V. CONCLUSIONS 

We have compared various gravitational-wave- 
extraction methods that arc nowadays very popular 
in numerical-relativity simulations: (i) the Abrahams- 
Price [65] technique based on the gauge-invariant 
Regge-Wheeler-Zerilli-Moncrief perturbation theory of 
a Schwarzschild space-time; (ii) the extraction method 
based on Weyl curvature scalars, notably the ^4 func- 
tion; (iii) some (variations of) quadrupole-type formulas. 
We have applied these methods to extract gravitational 
radiation from 3D numerical-relativity simulations of 
the very controlled system represented by a neutron star 
(with polytropic EOS), that is oscillating nonradially due 
to an initial pressure perturbation. The simulations have 
been performed via the Cactus-CCATIE-Carpet-Whisky 
general-relativistic nonlinear code. This code evolves 
the full set of Einstein equations in full generality in the 
three spatial dimensions. The accuracy of the waveforms 
extracted from the simulations, using the three methods 
recalled above, has been assessed (for small perturba- 
tions) via a comparison with waveforms (assumed to be 
exact) computed by means of the PerBACCo perturbative 
code. This code is designed to evolve, in the time 
domain, the Einstein equations linearized around a TOV 
background. It is 1+1-dimcnsional (i.e. one temporal 
and one spatial dimension) and adopts a constrained 
evolution scheme. This latter choice allows for the 
computation of very long and very accurate time series 
and similarly accurate waveforms. 

The initial pressure perturbation dp is given as an "ap- 
proximate" eigenfunction of the star, whose maximum is 
a fraction of the central TOV pressure p c . We focused 
only on £ = 2, m = 0, quadrupolar deformations, but 
we analyzed four values of the perturbation in order to 
cover the transition from the linear to the nonlinear os- 
cillatory regimes. We have first presented results of sim- 
ulations done using only the ID PerBACCo code to assess 
the accuracy of our exact waveforms. We have performed 



very long (about 1 s) and accurate simulations to extract 
both mode frequencies and damping times. We have an- 
alyzed finite-radius effects, finding that observers should 
be placed at extraction radius r > 200M in order to have 
amplitude errors below 1.6%. 

In doing 3D simulations in the perturbative regime, 
(10~ 3 ^ ma,x(8p/p c ) < 10~ 2 ), we have found that both 
metric and curvature wave-extraction techniques gener- 
ate waveforms that are consistent, both in amplitude and 
phasing, with the perturbative results. Each method, 
however, was found to have drawbacks. On one hand, the 
Zcrilli-Moncrief function presents an unphysical burst in 
the early part of the waveform; on the other hand, the 
ip4 scalar requires a polynomial correction to obtain the 
corresponding metric multipole. Our conclusion is that, 
in our setup, one needs both extraction methods to end 
up with accurate waveforms. 

For larger values of the initial perturbation amplitude, 
nonlinear effects in the 3D general-relativistic simula- 
tions are clearly present. The effective relative amplitude 
of the main modes of the extracted gravitational wave 
is smaller for larger amplitudes of the initial perturba- 
tion, because of mode couplings. The Fourier spectra 
of the rest-mass-density projections [see Eq. (78)] high- 
light that couplings between radial and quadrupolar fluid 
modes are present. Our study represents the first con- 
firmation, in fully general-relativistic simulations, of the 
results of Rcf . [96] , obtained via a perturbative approach. 

In addition, we have shown that the (non-gauge- 
invariant) generalizations of the standard Newtonian 
quadrupole formula that we have considered can be use- 
ful tools to obtain accurate estimates of the frequency of 
oscillation. By contrast, amplitudes are always signifi- 
cantly under/overestimated, consistently with precedent 
observations of Refs. [5, 85]. 

Finally, we discussed in detail some systematic errors 
that occur in the early part of the waveform extracted a la 
Abrahams-Price. These errors show up, in the early part 
of the Zerilli-Moncricf function, in the form of a burst of 
junk radiation whose amplitude grows linearly with the 
extraction radius. We have proposed some heuristic ex- 
planation of this fact and reproduced a similar behavior 
in low-accuracy perturbative simulations. Globally, our 
conclusion is that the extraction of the Zerilli-Moncrief 
function from a numerical-relativity simulation can be a 
delicate issue: small errors can conspire to give totally 
nonsensical results. Typically, these errors will show up 
as parts of the waveform whose amplitude grows with the 
observer's radius. We have also implemented the gener- 
alized wave-extraction approach based on the formalism 
of Refs. [6, 55, 56, 70], without any evident benefit. Note, 
however, that these kind of problems encountered with 
the Abrahams-Price wave-extraction procedure (as well 
as with its generalized version) seem to appear specif- 
ically in the presence of matter. In binary black-hole 
coalescence simulations curvature and metric waveforms 
seem to be fully consistent [32]. This last remark leads 
us to suggest that the Abrahams-Price wave-extraction 
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technique, a "standardized" and very basic procedure 
and infrastructure that has been developed long ago (and 
tested at the time) for specific applications to black-hole 
physics, should be rethought and reanalyzed when the 
Einstein equations are coupled to matter. For this rea- 
son, in the presence of matter, since systematic errors 
could be hard to detect and are present already in the 
simplest cases, we strongly encourage the community to 
make use of both wave-extraction techniques (curvature 
as well as metric perturbations) and to be always pre- 
pared to expect inaccuracies in the metric waveforms. 
In addition, concerning the many advantages related to 
extracting the metric waveforms directly from the space- 
time, we believe that it is also urgent and important for 
the community to have reliable implementations of the 
Abrahams-Price technique based on the Sarbach-Tiglio- 
Martel-Poisson [6, 55, 56, 70] formalism. 
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